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ABSTRACT 


Habitat evaluation constitutes an important and 
fundamental step in the management of wildlife 
populations and conservation policy planning. 
Geographic information system (GIS) and species 
presence data provide the means by which such 
evaluation can be done. Maximum Entropy (MaxEnt) 
is widely used in habitat suitability modeling due to its 
power of accuracy and additional descriptive properties. 
To survey snow leopard populations in Qomolangma 
(Mt. Everest) National Nature Reserve (QNNR), Xizang 
(Tibet), China, we pooled 127 pugmarks, 415 scrape 
marks, and 127 non-invasive identifications of the animal 
along line transects and recorded 87 occurrences 
through camera traps from 2014-2017. We adopted the 
MaxEnt model to generate a map highlighting the extent 
of suitable snow leopard habitat in QNNR. Results 
showed that the accuracy of the MaxEnt model was 
excellent (mean AUC=0.921). Precipitation in the driest 
quarter, ruggedness, elevation, maximum temperature 
of the warmest month, and annual mean temperature 
were the main environmental factors influencing habitat 
suitability for snow leopards, with contribution rates of 
20.0%, 14.4%, 13.3%, 8.7%, and 8.2% respectively. 
The suitable habitat area extended for 7 001.93 km“, 
representing 22.72% of the whole reserve. The regions 
bordering Nepal were the main suitable snow leopard 
habitats and consisted of three separate habitat patches. 
Our findings revealed that precipitation, temperature 
conditions, ruggedness, and elevations of around 4 
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000 m a.s.l. influenced snow leopard preferences at 
the landscape level in QNNR. We advocate further 
research and cooperation with Nepal to evaluate 
habitat connectivity and to explore possible proxies of 
population isolation among these patches. Furthermore, 
evaluation of subdivisions within the protection zones of 
QNNR is necessary to improve conservation strategies 
and enhance protection. 


Keywords: Qomolangma National Nature Reserve; 
Snow leopard; MaxEnt; Habitat suitability assessment; 
Tibet 


INTRODUCTION 


Wildlife habitat is defined as the surrounding environment where 
wild animals can accomplish their life cycle (Cody, 1987; Jiang et 
al., 2012). Habitats supply resources for population persistence, 
representing a determining factor for survival and successful 
reproduction (Wang & Chen, 2004; Yang et al., 2000). The 
habitat suitability index (HSI) is a measure of the ability for a 
habitat to sustain a species and is an important indicator for the 
quality of a given habitat (Jin et al., 2008; Lu et al., 2012; Song 
et al., 2014). Such evaluation constitutes one of the first steps in 
wildlife protection and management, offering a scientific rationale 
for the improvement of conservation policies (Liu et al., 2013). 
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The snow leopard (Panthera uncia) is a feline species 
distributed over 12 countries in Central Asia. It is estimated 
that China contains approximately 60% of the potential habitat 
available to snow leopards, who are reported to reside in 
the western provinces of Xinjiang, Xizang (Tibet), Qinghai, 
Gansu, Sichuan, and Inner Mongolia (McCarthy & Chapron, 
2003; Riordan & Shi, 2010, 2016). The species is currently 
classified as Vulnerable by the IUCN (International Union for 
Conservation of Nature and Natural Resources, 2017) and 
listed as a Class | protected animal by the China Key List 
(Riordan & Shi, 2010; Wang, 1998). Snow leopard-oriented 
research spans a diverse range of areas, including abundance 
and density (Alexander et al., 2015; Jackson et al., 2006), 
home range (Johansson et al., 2016), diet (Wang et al., 2014), 
behavior (Li et al., 2013), genetic diversity (Janecka et al., 
2017), climate change impact (Aryal et al., 2014d, 2016), 
human-snow leopard conflict (Aryal et al., 2014b; Chen et al., 
2016), and translocation of prey species (Aryal et al., 2013). 

Previous studies on snow leopard habitat use indicate two 
broad components determining such selection. First, snow 
leopard occurrence is predicted by several abiotic factors such 
as terrain (Slope and aspect) (Sharma et al., 2015; Wolf & Ale, 
2009), elevation (Alexander et al., 2016a; Aryal et al., 2014c), 
snow cover (Aryal et al., 2014c), distance from rivers (Aryal et 
al., 2014c), and annual mean temperature (Li et al., 2013). The 
second component is represented by biotic factors such as prey 
availability (Alexander et al., 2016b; Aryal et al., 2014a; Xu 8 
Luo, 2010) and human activity (Wolf 8 Ale, 2009). During the 
1990s, slightly more than 100 snow leopards were estimated to 
inhabit Qomolangma National Nature Reserve (QNNR). Based 
on three brief surveys, Jackson et al. (1994) estimated that 
“good” habitat totaled approximately 8 000 km?, or about 25% 
of QNNR's area. Furthermore, in May—June 2014 and May 
2015, a preliminary snow leopard presence survey (Chen et 
al., 2017) and a human-carnivore conflict survey (Chen et 
al., 2016) were conducted in QNNR, respectively. A more 
comprehensive and interdisciplinary study was also conducted 
to provide an evidential basis for the formulation of effective 
conservation policies and programs (Chen et al., 2017). Based 
on the above studies, habitat suitability assessment was 
deemed necessary for a systematic survey of the snow leopard 
population in QNNR. Several studies have been conducted 
to estimate the extent of suitable habitat covering the whole 
snow leopard range (McCarthy et al., 2016; Hunter & Jackson, 
1997; Fox, 1994). However, habitat suitability assessment at 
the regional level has not yet been reported, resulting in a gap 
between research and local conservation actions. 

Recently, with the develooment of 3S (GIS, RS, GPS) 
techniques, multiple models have been used to assess suitable 
habitat distribution, including mechanism models (Ouyang et 
al., 2001; Xu et al., 2006), regression models (Schadt et 
al., 2002), and ecological niche models (Su et al., 2015; 
Xu & Luo, 2010; Wang et al., 2008; Phillips et al., 2006). 
The maximum entropy (MaxEnt) model is an ecological niche 
model, which has been increasingly used to assess wildlife 
habitat distribution (Clements et al., 2012; Wilting et al., 2010). 
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MaxEnt models were originally formulated to estimate the 
presence density of a target species across a landscape 
(Phillips et al., 2006), but have since been applied to model 
species distribution and environmental niches, relying on 
“oresence-only” data and environmental predictors, thereby 
avoiding bias and leading to accurate results (Merow et al., 
2013; Pearson et al., 2007; Phillips et al., 2006). MaxEnt 
models are reliable, straightforward, and allow data to be easily 
obtained (Merow & Silander, 2014; Phillips & Dudik, 2008; 
Radosavljevic & Anderson, 2014). These models are suitable 
for the prediction of species habitat distribution at the whole 
reserve level, and an appropriate model for habitat suitability 
analysis (Haegeman & Etienne, 2010; Stachura-Skierczynska 
et al., 2009; Xing & Hao, 2011). 

Understanding the distribution of snow leopard habitat in 
QNNAR is important for the improvement of research outcomes 
and conservation plans. Thus, our study aimed to: (1) assess 
the potential habitat distribution for snow leopard in QNNR, 
along with biotic and abiotic factors of influence using the 
MaxEnt model; and (2) identify critical areas for snow leopard 
conservation under the existing functional zones of QNNR. 


MATERIALS AND METHODS 


Study area 

Located in the southwest Xizang (Tibet) Autonomous Region, 
China, QNNR (N27°48’-N29°19’, E84°27’—E88°23’) was 
established in 1989 to protect wildlife and ecosystems along 
the border of China and Nepal. The reserve covers an area 
of 33 814 km”, centering on the world’s highest peak, Mt. 
Everest. Altitude ranges from 1 440 m to 8 844 m a.s.l.. The 
average annual temperature is 2.1 °C and total annual rainfall 
is 270.5 mm. Some 81 species of mammals, 342 birds, 29 
amphibians and reptiles, and 8 fish inhabit the reserve (QNNR 
Administration, unpublished). Large mammals include the 
snow leopard (Panthera uncia), wolf (Canis lupus), lynx (Lynx 
lynx), brown bear (Ursus arctos), leopard (Panthera pardus), 
blue sheep (Pseudois nayaur), wild ass (Equus kiang), 
Tibetan gazelle (Procapra piticaudata), and Himalayan tahr 
(Hemitragus jemlahicus) (Chen et al., 2017). QNNR consists 
mainly of two broad ecosystems: that is, semi-humid mountain 
forest and semi-arid shrub along the southern and northern 
parts, respectively. Within these, vegetation varies across 
different sub-ecosystems. As altitude increases, 10 distinct 
sub-ecosystems can be observed progressively, including 
subtropical evergreen broad-leaved forest, subtropical 
semi-green broad-leaved forest, subtropical evergreen 
coniferous forest, warm temperature green coniferous forest, 
warm temperature sclerophyllous evergreen broad-leaved 
forest, subalpine temperature zone evergreen coniferous 
forest, subalpine temperature zone broad-leaved deciduous 
forest, alpine sub-frigid zone shrubland and grassland, alpine 
sub-frigid zone periglacial and alpine snow zone (QNNR 
Administration, unpublished). The perpendicular width of each 
vertical ecological system ranges from hundreds to thousands 
of meters, and each is very sensitive to external environment 
change (QNNR Administration, unpublished). For the different 


functional zones in QNNR, any disturbance is restricted in the 
core zone, scientific research can be conducted in the buffer 
area, and some human activity can occur in the experimental 
zone. 


Model selection 

Information science is the basis of MaxEnt models, which are 
widely used in many academic disciplines, as proposed by 
Jaynes (1957a, 1957b) and more recent studies (Jaynes & 
Bretthorst, 2003; Xing € Hao, 2011). The MaxEnt model is 
based upon ecological niche theory, relying on information from 
species “presence” data to explore the possible distribution of 
a target species within a study area. In 2004, MaxEnt software 
was developed to assess and evaluate suitable habitat for 


target species with good predictive power (Phillips et al., 2006). 


Occupancy models and generalized linear models (GLM) are 
also used to explore the relationship between species habitat 
selection and environmental factors (Alexander et al., 2015; 
MacKenzie et al., 2006; Wang et al., 2018). However, GLM 
and occupancy models exhibit similar limitations, such as the 
need for presence/absence or detection/non-detection data, 
which can introduce difficulties when an expansive study area 
is needed to detect rare species like snow leopards (Alexander 
et al., 2016a; MacKenzie et al., 2006). The MaxEnt model 
provides a self-checking function with automatically generated 
receiver operating characteristic (ROC) curves. Furthermore, 
to understand the regional level habitat distribution of snow 
leopards in QNNR, presence records are more likely to be 
accessible, and thus the MaxEnt model has an advantage 


over other types of models (Su et al., 2015; Liu et al., 2013). 


Therefore, the MaxEnt model was selected in this study to 
explore the potential suitable habitat of snow leopards in 
QNNR. 


Snow leopard “presence” data 

In 2014, a preliminary survey was conducted in four areas with 
33 camera traps around the villages of Zhalong (Jilong County), 
Dacang (Dingjie County), Riwu and Qudang (Dingri County) 
(Chen et al., 2017). Other study areas were later selected 
to estimate snow leopard abundance and density with high 
intensity. Our study area was divided into several 4 kmx4 km 
grids. We selected sampling transects and camera locations 
according to our knowledge of the ecological requirements 
of snow leopards and the experience of our local guide to 
maximize the detectability of animals in each grid. In 2015, 
following the preliminary survey, we selected Zhalong (Jilong 
County) and Qudang (Dingri County) as study areas. In total, 
83 and 23 camera traps (Ltl Acorn) were systematically set 
up over 400 km? in the Zhalong study area and 208 km? in 
the Qudang study area between October and November 2015, 
respectively. In summer 2016 (April to May), 112 camera traps 


(Ltl Acorn) were set up over 480 km? in the Zhalong study area. 


In winter 2016, an additional 68 camera trap sites were set up 
over 790 km? in another part of the Jilong study area. In winter 
2017, a snow leopard survey was conducted over 750 km? in 
the Zhaxizong study area (Dingri County). We recorded snow 
leopard signs (scats, tracks, scent marks, scrapes, bedding, 


and hair) along the most likely routes (ridgelines, hillsides, 
prey resource-rich places, and higher rugged places). Camera 
stations were placed with a minimum spacing of 1 km to 
maximize the number of individuals caught and adequately 
recapture individuals at different camera traps (Alexander et 
al., 2015). We placed two cameras in each site to capture both 
sides of the animal (Jackson et al., 2006). Within each grid, two 
or more camera stations were established. The Zhalong study 
area (Jilong County) was expanded from 400 km* in 2015 to 
790 km? in 2016. The transects used in 2015 were fixed for 
later surveys, but camera trap sites were not (Figure 1). 

For non-invasive genetics, we collected samples that 
belonged to snow leopards along transects to set up camera 
traps. The selection of feces was based upon shape, hair 
content, and terrain type or was found in association with 
other marking signs particular to this species. We collected 
scat samples using silica gel drying agent (Boitani & Powell, 
2012; Wasser et al., 1997). In the field, tubes containing 
silica desiccant were stored away from light in a dry and cool 
environment and were immediately stored at —20 °C after 
arrival at the laboratory. DNA from scat samples was extracted 
using a QlAamp DNA Stool Mini Kit (Qiagen, USA) following 
a modified extraction protocol under a dedicated UV-light 
laminar flow cabinet, physically separated by the pre-PCR 
area. Each batch of extracted samples was processed along 
with a negative control to account for possible contamination, 
with results screened on 1% agarose gel. Species 
identification was initially conducted through amplification of 
a 110 base pairs-portion of the cyt b (cyt b) gene using 
carnivore-specific primers, as described in Farrell et al. (2000) 
(forward: 5'-TATTCTTTATCTGCCTATACATACACG-3”; reverse: 
5’-AAACTGCAGCCCCTCAGAATGATATTTGTCCTCA-3’). PCR 
analysis was conducted under a dedicated UV-sterilized 
laminar flow hood using Premix Tag (Ex Tag Version 2.0; Takara 
Biotechnology) to a final volume of 20 uL containing 0.2 mmol/L 
of each dNTP, 2 mmol/L of Mg**, 0.4 umol/L of each primer, 7.4 
uL of DNAase/RNAase-free ultrapure water (Tiangen Biotech, 
Beijing), and 1 uL of extracted DNA. PCR conditions consisted 
of an initial denaturation step of 10 min at 94 °C, 35 cycles at 94 
°C for 30 s, 57 °C for 45 s, and 72 °C for 45 s, followed by a final 
elongation step at 72 °C for 10 min. All amplifications included 
One positive snow leopard reference sample and a blank 
control to account for contamination. Results were screened 
by 1% agarose gel electrophoresis. Successful PCR products 
were sequenced through an Applied Biosystems ABI 3730XL 
system by SinoGenoMax Limited Company (Beijing, China). 
Readings were assembled in Geneious 10.1 (Biomatters Ltd.) 
and matched against the NCBI database using BLAST to 
obtain the percentage of similarity. Only completely assembled 
sequences were considered for the dataset. Samples yielding 
incomplete assemblies, along with unsuccessful PCRs for cyt 
b and non-carnivore species identifications, were re-screened 
targeting a 126 base pairs-portion of the ATP6 gene using the 
carnivore-specific primers described in Chaves et al. (2012) 
(DF3: 5’-AACGAAAATCTATTCGCCTCT-3’; DR1: 5/-CCA 
GTATTTGTTTTGATGTTAGTTG-3’). The primers described in 
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Farrell et al. (2000) were used to target prey species DNA 
(Chaves et al., 2012). PCR volumes were identical to the 
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Figure 1 Location of camera traps for snow leopard surveys in Qomolangma National Nature Reserve, Xizang (Tibet) 


A: Qudang Study Area, Dingri County, winter 2015. B: Zhalong Study Area, Jilong County, winter 2015. C: Zhalong Study Area, Jilong County, summer 2016. D: 


Zhalong Study Area, Jilong County, winter 2016. 


Feces, scent marks, and killings could be attributed to 
carnivores other than snow leopards; however, pugmarks and 
scrapes are easily recognized. Therefore, pugmarks, scrapes, 
and feces identified as belonging to snow leopards were 
selected to create a suitable habitat distribution simulation in 
QNNR (Jackson & Hunter, 1996; Janecka et al., 2008). We 
used 1-km* grid cells for the MaxEnt niche-based modeling 
(Phillips et al., 2006). One snow leopard presence per cell 
was used and repeated location data were removed. A total 
of 222 GPS snow leopard presence points were included for 
further analysis. However, due to the higher number of field 
studies in the Jilong area, records were heavily geographically 
unbalanced, with the Jilong area containing more than 50% of 
all records. We thus further reduced the number of records 
in the Jilong area by randomly selecting records to produce 
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the same survey effort as outside the Jilong area. As only 59 
records were included for analysis outside the Jilong area (east 
QNNR), we only included 59 records from the Jilong area (west 
QNNR). 


Environment variables 

We tested the relationship of snow leopard presence data 
with possible influencing factors. We identified two main 
categories related to suitable habitat distribution, namely 
abiotic factors (elevation, slope, aspect, ruggedness, land 
use type, and distance factors) and bioclimatic and biotic 
factors (prey species, human related factors) (Alexander et al., 
2016a, 2016b; Sharma et al., 2015). We selected elevation, 
slope, aspect, ruggedness, land use type, and bioclimatic 
factors for landscape level analysis. We used snow leopard 
presence points and extracted 19 bioclimatic variables of our 


study area from WorldClim (www.worldclim.org) data (Table 1). 


In ArcGIS, elevation, aspect, slope, and ruggedness were 
determined using a digital elevation model (DEM) layer, with 
all variables clipped to our study areas. Ruggedness was 
extracted based on the method of Riley et al. (1999). The 
Topographic Ruggedness Index (TRI) was used to express the 
elevation difference between adjacent cells (all eight first-order 
neighbors within a quadratic grid) of a 90-m-resolution DEM 
on a scale from 1 (level) to 7 (extremely rugged). Land 
use and land cover were extracted from “GLOBCOVER 2009” 


Table 1 Variables used for modeling 


SN Variable Description 
1 BIO1 ( Continuous 
2 BlO2 (Mean diurnal range) Continuous 
3 BIOS (Isothermality) Continuous 
4 BlO4 (Temperature seasonality) Continuous 
5 BIO5 (Max temperature of warmest month) Continuous 
6 BIO6 (Min temperature of coldest month) Continuous 
7 BIO7 (Temperature annual range) Continuous 
8 BlO8 (Mean temperature of wettest quarter Continuous 
9 BIO9 (Mean temperature of driest quarter) Continuous 
10  BIO10 (Mean temperature of warmest quarter) Continuous 
11  —BlO11 (Mean temperature of coldest quarter) Continuous 
12 BlO12 (Annual precipitation) Continuous 
13 BlO13 (Precipitation of wettest month) Continuous 
14 BlO14 (Precipitation of driest month) Continuous 
15 BlO15 (Precipitation seasonality) SE 
(Coefficient of variation) 
16 BIO16 (Precipitation of wettest quarter) Continuous 
17 BIO17 (Precipitation of driest quarter) Continuous 
18 BIO18 (Precipitation of warmest quarter) Continuous 
19 BIO19 (Precipitation of coldest quarter) Continuous 


obtained from the University of Louvain (UCLouvain) and 
European Space Agency (ESA). Distance factors (distances 
to roads, settlements, and rivers) and biotic factors were 
excluded due to the small-scale snow leopard and prey species 
presence data. We then extracted the values of each variable 
corresponding to the presence locations of snow leopards to 
perform correlation analyses (Supplementary Table S1). After 
removing 10 highly correlated variables (>0.85), we used the 
remaining 14 variables for final analysis. 


Source 


Annual mean temperature) 


http://www.worldclim.org/bioclim 
(Hijmans et al., 2005) (1-19) 


Categorical (Irrigated croplands, Rainfed 


croplands, Mosaic Croplands/Vegetation, Mosaic 


Vegetation/Croplands, Closed to open broadleaved 


evergreen or semi-deciduous forest, Closed 


broadleaved deciduous forest, Closed needle-leaved 
evergreen forest, Closed to open mixed 


20 Land cover 


Forest-Shrubland/Grassland, Mosaic Grassland, 
Forest-Shrubland, Closed to open shrubland, 


broadleaved and needle-leaved forest, Mosaic 


GLOBCOVER 2009 
form University of 
Louvain (UCLouvain) 
and European Space 
Agency (ESA) (20) 


Closed to open grassland, Sparse vegetation, 
Close to open vegetation regularly 
flooded, Bare areas, Water 


bodies, Permanent snow and ice) 
Categorical (Flat; North: 0°—22.5°; 
Northeast: 22.5°-67.5°; East: 67.5°-112.5°; Southeast: 


21 Aspect 


112.5°-157.5°: South: 157.5°—202.5°:; Southwest: 


202.5°-247.5°:; West: 247.5°-—292.5°: Northwest: 
292.5°-337.5°; North: 337.5°—360°) 
Categorical (0-4.375, 4.375-9.310, 9.310-14.216, 


22 Slope 


GMTED2010 (21-24) 


14.216-19.223, 19.223-24.395, 24.395-—29.823, 


29.823-35.773, 35.773-43.812, 43.812—72.092) 


23 Elevation 


24 Ruggedness Continuous 


SN: Serial number. 


Continuous (m) 
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Simulation procedure 

Snow leopard presence data and selected variables were 
adapted to the format required for MaxEnt software (v 3.3.3k) 
(Phillips et al., 2006). We selected 75% of snow leopard 
presence data to build the model, with the remaining 25% used 


for model verification. We included 10 replicates in our analysis. 


We used a jackknife estimator to detect the importance of each 
variable. Sensitivity analysis was done for each variable with 
logistic output format. Results of the MaxEnt model were 
verified by ROC values: that is, rejected with a ROC value 
0.5-0.6; poor with 0.6-0.7; normal with 0.7-0.8; good with 
0.8-0.9; and excellent with 0.9-1.0. According to the expert 
experience method (Swets, 1988), the output results were used 
in ArcGIS 10.2 to reclassify the suitable snow leopard habitat 
distribution map. Using the MaxEnt model to evaluate habitat 
suitability in QNNR, the ASCII format file was imported into 
ArcGIS 10.2 for transformation into floating raster data. The 
floating raster was reclassified into low grade habitat (0-0.14), 
moderately suitable habitat (0.14—0.42), and highly suitable 
habitat (0.42-1). Focal statistics in ArcGIS were used to 
smooth the raster map to obtain the suitable snow leopard 
habitat distribution map in QNNR. 


RESULTS 


Snow leopard presence data records 

In 2014, we recorded five snow leopard pugmark sites and 65 
scrape sites, and 17 camera trap sites captured snow leopard 
photographs. In winter 2015, we recorded 38 pugmark sites, 
111 scrape sites, and 27 camera trap sites with photographs. In 
summer 2016, we recorded 35 pugmark sites, 131 scrape sites, 
and 44 camera trap sites with photographs. In winter 2016, we 
recorded 39 pugmark sites, 73 scrape sites, and 32 camera 
trap sites with photographs. In winter 2017, we recorded 10 
pugmark sites and 35 scrape sites. 

We collected a total of 52, 84, 135, and 72 fecal samples in 
summer 2014, winter 2015, summer 2016, and winter 2016, 
respectively. In 2014, we positively identified 29 samples 
through cyt b, 21 of which belonged to snow leopards. The 23 
unsuccessful samples, including two incomplete assemblies, 
produced 15 snow leopard identifications through the ATP6 
gene. In total, we identified 84% of samples (44), with 
snow leopard identifications accounting for 69% of the total 
(36). Of the 86 samples collected in winter 2015, 43 were 
positively identified as a carnivore species using cyt b, of 
which 14 samples belonged to snow leopard. The remaining 
43 were screened using the ATP6 gene; however, out of 
the positive PCRs, identification was only possible for one 
sample, which belonged to a snow leopard. In total, for 
samples collected in 2015, we identified 44 samples (51%), 
15 of which were snow leopards (17% of total). Species 
identification of the 135 samples collected in spring 2016 
yielded 71 assembled predator species identifications, 23 
of which were snow leopards. Out of the 64 samples 
selected for a second amplification, 55 failed to amplify, five 
belonged to blue sheep (Pseudois nayaur), one to yak (Bos 
grunniens), and three were partial assemblies belonging to 
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a carnivore species (either forward or reverse strand). We 
then amplified 64 samples for the ATP6 gene, yielding 31 
positive carnivore identifications, with 26 belonging to snow 
leopards. In summary, 102 samples (75.5%) were successfully 
identified, with snow leopards representing 36.2% of the total 
(49 identifications). For winter 2016, we identified 24 carnivores 
out of 72 samples through cyt b, 11 of which were identified as 
snow leopard. The remaining 48, amplified through the ATP6 
gene, yielded 19 carnivores, with snow leopards identified 16 
times. In total, we identified 43 carnivores in 43 occasions 
(60% of total), and snow leopards represented 37.5% of total 
fecal samples. Pooling all samples together over the four field 
seasons, we collected 345 fecal samples and identified 233 
(66%) as belonging to a predator species. Snow leopards were 
identified 127 times, representing 36.8% of the total collected 
scats and 54.5% of the successful identifications. Results are 
summarized in Table 2. 


Table 2 Identification results from molecular scatology 


Field Season Scats Cytb+ ATP6+ SLs Failed 


Summer 2014 52 29 15 36 8 
Winter 2015 84 43 1 15 42 
Summer 2016 135 71 31 49 33 
Winter 2016 72 24 19 27 29 
Total 343 167 66 127 112 


Cyt b+ and ATP6+ indicate success of identification using the two 
markers. SLs: Snow leopards. 


In total, 127 pugmark sites, 415 scrape sites, and 127 molecular 
identifications of snow leopard were recorded. In addition, 120 
camera trap sites captured snow leopard images (Figure 2). 


MaxEnt prediction evaluation 

The ROC results (Figure 3) showed a mean AUC value of 
0.921, indicating that the predictions obtained from the MaxEnt 
model were excellent. 


Suitable snow leopard habitat distribution with 
environment factors 

According to the suitable habitat distribution map, we 
determined that suitable habitat was located mainly along 
mountain areas near the Nepalese border. We identified three 
distinct separated patches: (1) Zhalong and Gongdan areas 
in Jilong county; (2) Chabuling area in Nielamu county and 
Rongxia area in Dingri county; and, (8) Qudang area and 
perimeter zone in Dingri county. 

The jackknife estimator results (Figure 4; Table 3) showed 
that precipitation in the driest quarter (BIO17; 20.0%), 
ruggedness (14.4%), elevation (13.3%), maximum temperature 
of the warmest month (BIO5; 8.7%), and annual mean 
temperature (BIO1; 8.2%) were the main factors contributing 
to snow leopard habitat selection. The importance rates in the 
MaxEnt model prediction indicated that ruggedness (34.4%), 
mean diurnal range (BIO2; 16.9%), maximum temperature 
of the warmest month (BIO5; 11.3%), annual precipitation 
(BIO12; 8.8%), and precipitation in the driest quarter (BlO17; 


7.9%) were the five main factors affecting snow leopard habitat 
preferences (Table 3; Table 4). 

Sensitivity analysis determined the influence of each factor 
on snow leopard habitat distribution (Figure 5). The altitude in 
areas suitable for snow leopard habitat was between 3 900 m 


and 5 000 m a.s.l.. The ruggedness range most suitable 
for snow leopards was between 1 000 m and 1 300 m a.s.l. 
where the terrain was extremely rugged. Consequently, snow 
leopards preferred to use more rugged terrain at elevations of 
around 4 000 m a.s.l.. 








Figure 2 Snow leopard presence data collected in QNNR 
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Figure 3 ROC verification of distribution of suitable snow leopard habitat in QNNR 
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Jackknife of regularized training gain for snow leopard 


Aspect 
Bio_1 
Bio_12 
Bio_13 
Bio_14 
Bio_17 
Bio_2 
Bio_3 
Bio_5 
Bio_6 


Environmental variable 


Elevation 
Landcover 
Ruggedness 


Slope 


0.0 0.2 0.4 0.6 


Without variable * 
With only variable ™ 
With all variables ® 


1.0 1.2 1.4 


Regularized training gain 


Figure 4 Jackknife test of environmental variables in training data by MaxEnt 


Distribution of suitable snow leopard habitat in QNNR 


From the reclassified map, the area of each habitat class 
(low grade habitat, moderately suitable habitat, and highly 
suitable habitat) was calculated. The area of highly suitable 
habitat in QNNR was 1 756.55 km*, moderately suitable habitat 
was 5 245.38 km*, and low-grade habitat was 23 814.05 
km?. Thus, the area of relatively good snow leopard habitat 
totaled 7 001.93 km? in QNNR, accounting for 22.72% of the 
overall area of QNNR (Figure 6). According to the existing 
functional zones (core, experimental, and buffer zones) of 
QNNR (Figure 7), only 23.24% (1 627.52 km*/7 001.93 km?) 
of good habitat lies in the core zone, 36.06% (2 525.22 
km?/7 001.93 km?) lies in the buffer zone, and 40.52% 
(2 837.40 km?/7 001.93 km?) lies in the experimental zone. 


Table 3 Contribution and permutation importance values of 
environmental variables 


Environmental variables Contribution (%) Permutation importance (%) 





Aspect 1.3141 1.38 
BIO1 8.1786 0.4244 
BlO12 6.8994 8.8319 
BIO13 4.214 2.6905 
BlO14 0.7669 2.2132 
BlO17 19.9978 7.8506 
BlO2 6.3171 16.8578 
BIO3 4.4267 3.8942 
BIO5 8.7044 11.2613 
BIO6 4.1076 0.508 
Elevation 13.3151 1.6187 
Land cover 2.3265 3.9645 
Ruggedness 14.3777 34.3892 
Slope 5.154 4.1157 
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DISCUSSION 


Environmental factors influencing snow leopard habitat 
selection 

Our results showed that precipitation in the driest quarter (BIO17), 
ruggedness, elevation, maximum temperature of the warmest 
month (BIO5), and annual mean temperature (BIO1) were the five 
main factors influencing snow leopard habitat suitability in QNNR. 
Wolf & Ale (2009) conducted research in the Sagarmatha National 
Park (area of 1 148 km?) of Nepal and reported that terrain and 
human activity were the main factors determining snow leopard 
spatial distribution, whilst prey species had a moderate effect. 
In Nepals Annapurna Conservation Area, Aryal et al. (2014c) 
indicated that cliffs, grassland, and shrubland at high elevations (3 
000-5 000 m a.s.l.) were preferred habitats of snow leopards in 
the study area (about 2 025 km?). A study conducted in Mongolia 
reported that distribution of prey resources and rugged terrain 
largely explained changes in snow leopard habitat use (McCarthy et 
al., 2005). In India, investigation in intensively grazed areas showed 
that snow leopard habitat-use mainly relied on wild prey species 
density (Sharma et al., 2015). In China, very little research has been 
conducted on snow leopard habitat use. A winter habitat-use survey 
of snow leopards in Tomur National Nature Reserve highlighted 
that prey resources and principal geographic features (ruggedness, 
bases of cliffs, and stream beds) were the main factors influencing 
snow leopard habitat use within the 3 000 km? study area (Xu et al., 
2012). In Sanjiangyuan National Nature Reserve, a landscape level 
analysis of snow leopard habitat using the MaxEnt model indicted 
that annual average temperature and ruggedness were the two 
main factors influencing habitat selection (Li et al., 2013). From 
the above studies, we can conclude that the determining factors of 
snow leopard habitat selection may differ in different areas. Thus, 
local surveys on snow leopard habitat selection are critical to adapt 
conservation needs according to the local context. 


Table 4 Cumulative and logistic thresholds and corresponding omission rates used for modeling 


Fractional Trainin 
Cumulative threshold Logistic threshold Description g 


predicted area omission rate 


1.000 0.015 Fixed cumulative value 1 0.641 0.000 
5.000 0.056 Fixed cumulative value 5 0.382 0.014 
10.000 0.113 Fixed cumulative value 10 0.262 0.047 
3.590 0.042 Minimum training presence 0.464 0.000 
18.200 0.203 10 percentile training presence 0.173 0.099 
24.048 0.261 Equal training sensitivity and specificity 0.128 0.128 
30.946 0.330 Maximum training sensitivity plus specificity 0.096 0.133 
5.371 0.060 Balance training omission predicted area 0.376 0.006 
and threshold value 
Equate entropy of thresholded and 
13.893 0.159 A 0.210 0.086 
original distributions 
Aspect Bio_1 Bio_2 
1.0 1.0 1.0 
0.0 0.0 0.0 
-1 359.869 -175 148 102 167 
Bio_3 Bio_5 Bio_6 
1.0 1.0 1.0 
0.5 0.5 0.5 
0.0 0.0 0.0 
44 40 59 238 331 22 
Bio_12 Bio_13 Bio_14 
1.0 1.0 1.0 
0.5 a 0.5 TN 0.5 
0.0 0.0 0.0 
206 2587 71 893 0 13 
Bio_17 Elevation Landcover 
1.0 1.0 1.0 
0.5 0.5 0.5 | ene 
0.0 0.0 0.0 
3 73 1913 7860 11 220 
Ruggedness Slope 
1.0 1.0 
0.5 0.5 
0.0 0.0 
-0 2448.403 0 86.448 


Figure 5 Response curve of selected variables for snow leopard habitat suitability in ANNR 
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Figure 6 Distribution of suitable snow leopard habitat in QNNR 
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Figure 7 Distribution of suitable snow leopard habitat in different functional zones in QNNR 
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Little has been done to understand regional level snow 
leopard habitat suitability. Our results found that precipitation 
and temperature conditions at the regional level had a strong 
influence on suitable snow leopard habitat. Thus, it is likely that 
climate change will influence snow leopard habitat selection (Li 
et al., 2016). Among abiotic factors, elevation and ruggedness 
had a greater influence on habitat suitability. Our study area 
had the highest average elevation among the whole global 


range of snow leopards, ranging from 3 500 m to 5 500 m a.s.!.. 


Our results indicated that snow leopards preferred elevations of 
around 4 000 m a.s.!. near Mt. Everest. Like previous research, 
our findings confirmed that snow leopards favored highly to 
extremely rugged areas, as based on the highest permutation 
importance values (Li, 2013). Highly rugged areas usually 
contain large-sized rugged rocks, which may provide shelter 
for snow leopards, and thus represent an important feature for 
their survival. The response curves of land cover and slope 
showed that snow leopards preferred bare and relatively high 
slope areas. Aspect had little influence on snow leopard habitat 
suitability in QNNR. 

We selected the MaxEnt model as the most appropriate 
method to accomplish our research aims. However, we 
recognize that bias may still exist due to the relatively small 
study area and the inclusion of a limited number of habitat 
types, mainly due to resource availability and accessibility. This 
highlights the common practice of research groups to focus 
economic efforts toward survey areas where the probabilities 
of encountering a target species are higher. One fundamental 
assumption of the MaxEnt model is that the entire area of 
interest has been systematically or randomly sampled (Phillips 
et al., 2009; Royle et al., 2012). Our survey showed a strong 


sampling bias toward some regions or environmental features. 


To account for this factor, we selected 222 presence sites 
according to the pixel cells of environmental variables. Many 
occurrence records were available for this study, therefore we 
adopted spatial filtering and balancing of occurrence data to 
minimize omission errors (false negatives) and commission 
errors (false positives) (Kramer-Schadt et al., 2013). Spatial 
clumping of records was reduced in datasets for MaxEnt model 
calibration. Environmental variables were tested regarding their 
autocorrelation relationship. The manipulation of background 
datasets was considered inaccurate, as it might increase the 
risk of omission errors for species like snow leopards with a 
generalist response in many predictors. 

Based on the collected data, we primarily selected non-biotic 
factors for further analysis. However, the limited number of 
examined factors may be inappropriate (at least to some extent) 
for our landscape level analysis in QNNR. An increase in 
the number of evaluated factors is preferable for landscape 
level investigation, as advocated by Robinson & Weckworth 
(2016). More meaningful variables (e.g., human-related factors) 
influencing snow leopard habitat distribution should be included 
in future research. Understanding snow leopard habitat in 
QNNR is important, thus this study provides a basis for more 
in-depth analysis on population densities and represents a 


starting point to implement wider conservation-oriented studies. 


Suitable habitat distribution 

The MaxEnt model was used to assess habitat suitability 
in QNNR and, according to the AUC values, results were 
excellent. We showed that suitable snow leopard habitat 
in QNNR was mainly located along the border with Nepal, 
with three distinct habitat patches detected within the nature 
reserve (Figure 6). This habitat separation might be due to the 
presence of very high mountains, including Mt. Everest and the 
Xixiabangma area (Figure 6). According to the geographical 
features and local knowledge (personal communication), snow 
leopards are reported to exist along the southern border of the 
Xixiabangma area and in the northern area of Mt. Everest (near 
base camp). Using the HIS method, Jackson et al. (1994) 
described preliminary results for suitable snow leopard habitat, 
showing such habitat to be mainly located in the western 
area of QNNR. Our results indicated that “good” snow leopard 
habitat area was about 1 000 km? smaller than the results given 
by Jackson et al. (1994). This might be a consequence of the 
increase in human population and activities like tourism (Chen 
et al., 2017), or because of the variables and methods used for 
habitat suitability estimation. 

Further study should focus on the three unconnected habitat 
patches to assess whether individual snow leopards are subject 
to isolation. Cross-boundary cooperative research with Nepal 
is also necessary. In QNNR, low-grade snow leopard habitat 
accounted for about 84% of all potential snow leopard habitat 
and included very flat areas with little shelter or usable rock 
cover. However, it is necessary to survey these areas to 
confirm the “absence” of snow leopards and understand natural 
prey population statuses. Low-grade snow leopard habitat is 
important for their survival because it constitutes a part of the 
ecosystem where snow leopards could exist. 

Good snow leopard habitat in QNNR mainly lies in 
the experimental zone, which is not optimal for protection. 
Therefore, assessing the efficiency of the functional zones in 
QNNR is needed in order to implement better conservation 
strategies and promote new protection policies. 

QNNR is located in the core area of the Himalayas and 
is an important area for global snow leopard conservation 
(McCarthy & Chapron, 2003). Since its establishment 30 years 
ago, no-grazing and controlled-grazing policies have been 
implemented in QNNR; however, the contradiction between 
protection and development is still prominent (Chen et al., 
2016). To improve protection of rare species and ecosystems in 
QNNR, it is imperative to strengthen management of the whole 
nature reserve, accounting for target species occurrence and 
human conflict. 
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ABSTRACT 


DNA damage in oocytes can cause infertility and birth 
defects. DNA double-strand breaks (DSBs) are highly 
deleterious and can substantially impair genome 
integrity. Homologous recombination (HR)-mediated 
DNA DSB repair plays dominant roles in safeguarding 
oocyte quantity and quality. However, little is known 
regarding the key players of the HR repair pathway 
in oocytes. Here, we identified oocyte-specific gene 
Ooep as a novel key component of the HR repair 
pathway in mouse oocytes. OOEP was required for 
efficient ataxia telangiectasia mutated (ATM) kinase 
activation and Rad51 recombinase (RAD51) focal 
accumulation at DNA DSBs. Ooep null oocytes were 
defective in DNA DSB repair and prone to apoptosis 
upon exogenous DNA damage insults. Moreover, Ooep 


null oocytes exhibited delayed meiotic maturation. 


Therefore, OOEP played roles in preserving oocyte 


quantity and quality by maintaining genome stability. 


Ooep expression decreased with the advance of 
maternal age, suggesting its involvement in maternal 


aging. 


Keywords: Ooep; Homologous recombination; DNA 
double-strand break repair; ATM; RAD51 


INTRODUCTION 


During the long period of meiotic arrest, oocytes are 
exposed to endogenous and exogenous genotoxic insults 
and tend to accumulate DNA damage. DNA damage can 
occur as single-stranded breaks (SSBs) or double-strand 
breaks (DSBs), the latter of which are highly deleterious 
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and can substantially impair genomic integrity and cell 
viability. Two major pathways are involved in DNA DSB 
repair, non-homologous end joining (NHEJ)- and homologous 
recombination (HR)-mediated repair. In oocytes, defection 
of the error-free HR pathway can cause a decline in oocyte 
quantity and quality in natural and pathological reproductive 
aging in both mice and humans (Day et al., 2015; Perry et 
al., 2013; Titus et al., 2013; Wang et al., 2014). Therefore, 
HR-based DNA DSB repair in oocytes plays a key role in 
regulating female reproductive performance. This repair has 
been intensively studied in somatic cells, with many players 
identified (Adamson et al., 2012); however, little information is 
available regarding the components of the HR repair pathway 
in oocytes. Identifying essential participants will shed light on 
understanding oocyte and ovarian aging. 

Oocyte-expressed protein (Ooep) is an oocyte-specific 
maternal-effect gene. The OOEP protein localizes at the 
subcortical region of an oocyte and interacts with four other 
maternal proteins, including MATER, FILIA, PADI6, and TLE6, 
to form a subcortical maternal complex (SCMC) in mouse and 
human oocytes (Li et al., 2008; Zhu et al., 2015). Depletion of 
OOEP or its interaction proteins can lead to embryo arrest at 
the 2-cell stage with uncharacterized molecular mechanisms in 
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mice (Li et al., 2008; Tong et al., 2000; Yu et al., 2014; 
Yurttas et al., 2008; Zheng & Dean, 2009). We recently 
utilized mouse embryonic stem cells as a model to explore the 
function of FILIA and uncovered its critical roles in safeguarding 
the genomic stability of pluripotent stem cells (Zhao et al., 
2015). This finding prompted our hypothesis that OOEP may 
be involved in regulating genome stability in oocytes. Thus, in 
this study, we carefully examined the function of OOEP in DNA 
DSB repair in oocytes. 


MATERIALS AND METHODS 


Animals 

The mouse line with targeted mutation of Ooep (background 
C57BL/B6xSv129) was kindly provided by Dr. Jurrien Dean 
from the National Institutes of Health, USA. These Ooepi”m 
female mice do not express the OOEP protein in oocytes (Li 
et al., 2008). Mice were maintained in specific pathogen-free 


conditions. All experimental procedures and animal care 
were performed according to the protocols approved by the 
Ethics Committee of the Kunming Institute of Zoology, Chinese 
Academy of Sciences. 


DNA damage treatment of oocytes and ovaries 


Germinal vesicle (GV)-stage oocytes were cultured according 
to standard procedures (Marangos & Carroll, 2012). Ovaries 
were collected from newborn mice at postnatal day 5 (P5) and 
cultured according to standard protocols (Gonfloni et al., 2009). 


Antibodies 


Rabbit anti-OOEP serum was raised against the 1-19 amino 
acids of the OOEP protein (Li et al., 2008), and antibody 
specificity was verified by utilizing Ooep null oocytes. Other 
primary and secondary antibodies were obtained commercially, 
with relevant information shown in Table 1. 


Table 1 Information on primary and secondary antibodies for immunofluorescence staining (IF) and immunoblotting (IB) 


Antibody Company Host Catalog number Dilution 
Primary 
1:200 (IF 
y-H2AX Cell signaling, USA Rabbit 9718S IF) 
1:800 (IB) 
1:200 (IF 
p-AT M(Ser1981) Novus biologicals, USA Mouse NB100-306 A 
1:800 (IB) 
1:200 (IF 
RAD51 Abnova, China Mouse H00005888-B01P a, 
1:500 (IB) 
DDX4 Abcam, USA Rabbit AB13840 1:200 (IF) 
Beta-ACTIN Zhongshan, China Mouse TAO9 1:5000 (IB) 
GFP Abcam, USA Chicken AB13970 1:2000 (IB) 
Secondary 
Alexa Fluor 488 Donkey anti-Rabbit IgG Invitrogen, USA A-21206 1:500 (IF) 
Alexa Fluor 555 Donkey anti Rabbit IgG Invitrogen, USA A-31572 1:500 (IF) 
Alexa Fluor 647 Goat anti Rabbit IgG Invitrogen, USA A-21244 1:500 (IF) 
Alexa Fluor 488 Goat anti Mouse IgG Invitrogen, USA A-11029 1:500 (IF) 
Alexa Fluor 555 Donkey anti Mouse IgG Invitrogen, USA A-31570 1:500 (IF) 
Alexa Fluor 488 Goat anti Chicken IgG Invitrogen, USA A-11039 1:500 (IF) 


Immunofluorescence 
quantification 
Immunofluorescence staining of oocytes and ovaries was 
performed according to standard procedures (Guo et al., 
2016; Xiong et al., 2008). DNA was labeled with 4’, 
6-diamidino-2-phenylindole (DAPI, Sigma, Germany). Images 
were analyzed using FV10-ASW 2.1 viewer software. The 
intensities of the foci were quantified using Image J software. 


staining and foci intensity 


Comet assay 

Alkaline comet assay was performed according to standard 
procedures (Berthelot-Ricou et al., 2011; Tice et al., 2000). 
Comets were analyzed using CASP comet assay analysis 
software (Andor Technology, UK). Three independent repeats 
were performed. 
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Detection of apoptosis by TUNEL staining 

Frozen ovarian sections were permeabilized with 0.1% 
Triton X-100 in phosphate-buffered saline (PBS), followed by 
incubation with the TUNEL reaction solution from the In-Situ 
Cell Death Detection Kit, Fluorescein (Roche Diagnostics, 
USA). All staining procedures were performed according to the 
manufacturer’s instructions. DNA was labeled with DAPI. 


Hematoxylin and eosin (HE) staining and primordial and 
primary follicle counting 

Ovarian sections were incubated with the HE reaction solution 
(BOSTER, USA) at room temperature (RT). All staining 
procedures were performed according to the manufacturer's 
instructions. Oocytes residing in primordial and primary follicles 
were counted, as described previously (Myers et al., 2004; 


Skaznik-Wikiel et al., 2007). 


Construction of OOEP-GFP expression plasmid 

The protein-coding region of Ooep was amplified by 
polymerase chain reaction (PCR) and inserted into 
pcDNA3.1/CT-GFP-TOPO (Invitrogen, USA). Integrity was 
confirmed by DNA sequencing. 


In vitro transcription and mRNA microinjection of GV 
oocytes 

Ooep-Gfp expressing plasmids were linearized with Scal (New 
England Biolabs, USA). The mRNA was synthesized using 
an in vitro transcription kit (mMessage mMachine T7 kit, 
Ambion, USA), and purified with a RNeasy MinElute Cleanup 
Kit (Qiagen, Germany). The mRNA was then dissolved in 
nuclease-free water and stored at —80 °C. We microinjected 
500 ng/uL of MRNA in injection buffer (10 mmol/L Tris-HCl (pH 
7.5) and 0.1 mmol/L EDTA) into the cytoplasm of Ooep-/- GV 
oocytes. 


Single-oocyte cDNA amplification and quantitative RT-PCR 
Briefly, cDNA was prepared from a single GV oocyte and 
amplified by 20 cycles of PCR according to published 
protocols (Tang et al., 2010). Glyceraldehyde 3-phosphate 
dehydrogenase (Gapdh) was used as the housekeeping 
control. Primers for amplification of Gapdh and Ooep included: 
Ooep forward: 5'-GTCATAGGCACAGACCAAGCG-3', Ooep 
reverse: 5’-GGCCGCCAT GT TCAAGAGAAT-3’; Gapdh forward: 
5’-TTGAGGTCAATGAAGGGGTC-3’, Gapdh reverse: 5’-TCG 
TCCCGTAGACAAAAT GG-3’. GraphPad Prism 5 software was 
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Figure 1 Ooep depletion causes DNA damage in mouse oocytes 





used for statistical analysis. 


Statistical analyses 

Quantitative analyses were based on at least three 
independent repeats and results were represented as 
means+SEM. Data were first tested by the homogeneity of 
variance test. For non-normal distribution, the data were 
subjected to nonparametric tests. Otherwise, the data were 
analyzed by t-tests. We considered P<0.05 as statistically 
significant. 


RESULTS 


Ooep participates in DNA double-strand break repair in 
mouse oocytes 

Upon DNA damage, histone H2AX is phosphorylated at 
Ser139 (y-H2AX) and recruited to the damaged sites to form 
visible foci under confocal microscopy (Rogakou et al., 1998). 
y-H»AX foci formation is generally considered as a marker of 
DNA damage. In the fully-grown GV oocytes, depletion of 
OOEP caused a significant increase in y-H2AX foci intensity 
compared to wild-type oocytes (Figure 1A), suggesting more 
endogenous DNA damage in Ooep*™™*™ oocytes. To validate 
this observation, we performed comet assay, an unambiguous 
method that measures the extent of DNA damage on a single 
cell basis (Berthelot-Ricou et al., 2011; Tice et al., 2000). 
The GV oocytes from Ooep'™” females displayed significantly 
longer comet tails than those from the wild-type counterparts 
(Figure 1B), confirming that Ooep!”''” oocytes contained more 
endogenous DNA damage. 





Olive moment 
NO 
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A: Immunostaining revealed higher levels of y-H2AX foci in Ooep!'”'” GV oocytes than in wild-type (WT) counterparts. Quantification of the y-H2AX foci intensity 


is shown in the lower panel. B: Comet assay confirmed that Ooep!”""” GV oocytes showed greater DNA damage. Data are presented as mean+SEM. Scale 


bars: 10 um, ***: P<0.001. 
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The accumulation of DNA damage in Ooept”m oocytes 
suggested inefficiency in DNA damage repair. To test this 
hypothesis, we treated wild-type and mutant GV oocytes 
with 50 ug/mL of etoposide, a topoisomerase Il inhibitor, 
to induce DNA DSBs (Marangos & Carroll, 2012; Nagy & 
Soutoglou, 2009), and compared the dynamics of y-HsAX 
resolution. Immediately after treatment, wild-type and mutant 
oocytes had comparable y-HsAX levels, as measured by 
immunostaining (Figure 2A) and immunoblotting analyses 
(Figure 2B). Following several hours of DNA repair recovery, 
y-H2AX was significantly resolved in the wild-type oocytes, 
reflecting efficient DNA damage repair. In sharp contrast, 
y-H2AX remained at a higher level in the Ooep’””™ oocytes 















than in the wild-type oocytes, indicating compromised DNA 
damage repair (Figure 2A, B). To further validate the 
role of OOEP in DNA damage repair, we performed a 
rescue experiment by micro-injecting Gfo (green fluorescent 
protein)-tagged Ooep mRNA into the Ooep’”™” oocytes. 
Consistently, the OOEP-complemented oocytes resolved 
y-H2AX more efficiently than the Ooeptmm oocytes after 
etoposide treatment and recovery (Figure 2C). These data 
together support the notion that maternal protein OOEP is 
necessary for efficient DNA DSB repair in oocytes. In line with 
this role, the protein expression of OOEP was slightly induced 
by etoposide treatment in GV oocytes (Figure 2D). 
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Figure 2 Ooep participates in DNA double-strand break repair in mouse oocytes 
A: Ooep"”" and WT oocytes were subjected to etoposide treatment and recovery. Immediately after etoposide treatment (no recovery), Ooep"””” and WT 
oocytes showed a similar level of y-H2AX. After 6-h recovery, Ooep"””” oocytes retained more y-H2AX foci than WT oocytes. B: Etoposide treatment induced 


tmtm oocytes after 4-h recovery. Each sample 


a similar increase in y-H2AX in both WT and Ooep!”'” oocytes. y-H2AX was better resolved in WT but not Ooep 
contained 80 oocytes. C: OOEP-GFP was re-expressed in Ooep!”'” oocytes (left panel). After 3-h recovery from the same etoposide treatment, OOEP-rescued 
oocytes had a lower level of y-H2AX than Ooep'” oocytes (right panel). D: Immunoblotting analysis showed that OOEP protein level were upregulated by 


etoposide treatment. Each sample contained 20 oocytes. Data are presented as meant SEM. Scale bars: 10 um. *: P<0.05. 


2015; Titus et al., 2013). ATM activation by DSBs initiates 
robust checkpoint signaling as well as DSB repair processes 
in somatic cells (Behrens et al., 2014; Jackson & Bartek, 2009) 
and germ cells (Di Giacomo et al., 2005; Oktay et al., 2015). 


Ooep may be involved in HR-mediated DNA DSB repair 


In oocytes, DNA DSB is predominantly repaired via the 
HR-mediated pathway, in which ataxia telangiectasia mutated 
(ATM) kinase plays a central orchestration role (Oktay et al., 
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In mouse GV oocytes, ATM can be activated by certain levels 
of etoposide-induced DNA DSBs. We therefore treated the 
oocytes with 50 ug/mL of etoposide for 3 h, as per previous 
study (Marangos & Carroll, 2012), and examined the effect of 
OOEP on AIM activation by phosphorylation at serine 1981 
residue (p-ATM) (Lavin 8 Kozlov, 2007). Co-immunostaining 
analysis showed that fewer p-ATM foci were formed and 
co-localized with y-H2AX foci in the Ooept”'T oocytes than in 
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the wild-type oocytes, despite y-H2AX level being comparable 
(Figure 3A). Immunoblotting analysis of p-ATM further validated 
this observation. Under normal conditions, very little p-ATM 
was detected in either wild-type or mutant oocytes. However, 
higher levels of p-ATM were induced in wild-type oocytes 
compared with mutant oocytes following etoposide treatment 
(Figure 3B). These results suggest that OOEP regulates ATM 
activation upon DNA DSBs. 
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Figure 3 OOEP is required for ATM activation and RAD51 recruitment to DNA damage sites 


A: Ooep'""™ and WT oocytes were subjected to the same etoposide treatment. Immunofluorescence staining showed that at 3-h post damage, Oovep!”'” 


oocytes contained a lower level of active ATM (p-ATM) compared with WT counterparts. B: Immunoblotting detection of p-ATM level under untreated and 


etoposide-treated conditions. In untreated conditions, rare p-ATM was detected in WT and Ooep""” oocytes. At 3 h after etoposide treatment, higher level of 


p-ATM were detected in WT but not Ooep'” oocytes. Each sample contained 80 oocytes. C: Immunofluorescence staining showed that at 3-h post etoposide 


treatment, Ooep'” oocytes had a lower level of RAD51 than WT controls. Data are presented as mean+SEM. Scale bars: 10 um. ***: P<0.001, *: P<0.05. 


RAD51 is an essential recombinase downstream of ATM 
in the HR pathway (Oktay et al., 2015, Perez et al., 2007). 
After DNA DSBs, the RAD51 protein is recruited and binds 
to resected single-strand DNA to form DNA-protein filaments 
(RAD51 foci), which facilitate the searching and pairing of DNA 
homologues during recombination (Jasin & Rothstein, 2013). 
Upon etoposide treatment, RAD51 foci were formed at damage 
sites labeled with y-H»AX in wild-type oocytes. However, fewer 
RAD51 foci were detected in Ooep'”m oocytes (Figure 3C). 


Thus, the absence of OOEP in oocytes impaired RAD51 focal 
accumulation at DNA damage sites. Taken together, these data 
demonstrate that OOEP may function in the HR-mediated DNA 
DSB repair pathway by regulating ATM activation and RAD51 
recruitment to DSB sites in mouse oocytes. 


Ooep protects oocytes from DNA DSB-induced apoptosis 
and meiosis delay 


Persistent DNA DSBs can evoke p63-mediated apoptosis of 
oocytes in primordial and primary follicles (Gonfloni et al., 
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2009; Suh et al., 2006). OOEP is necessary for HR-mediated 
DNA DSB repair in oocytes, suggesting that Ooep!”™ 
oocytes might be sensitive to exogenous or endogenous 
DNA damage insults and prone to apoptosis. To test 
this hypothesis, we compared the apoptosis susceptibility 
of oocytes to exogenous DNA damage insults between 


wild-type and Ooep!""™ females at the early postnatal stage. 


Five-day-old (P5) mouse ovaries, which contained mostly 
primordial follicles, were collected and cultured for 12 h with 
various doses of DNA cross-linking agent cisplatin (Gonfloni 
et al., 2009). At the concentration of 3 mg/L, cisplatin was 
able to induce mild apoptosis (~13.5%) in wild-type ovaries 
and was therefore utilized in the following studies. After in 
vitro treatment, ovarian sections from Ooep'™” and wild-type 
females were counterstained with DDX4 (a germ cell specific 
marker) and TUNEL. The percentages of apoptotic oocytes 
(DDX4*TUNEL*) were significantly higher in Ooep'"7 infants 
than in wild-type counterparts (Figure 4A), suggesting that 
Ooep!t"imM oocytes were more sensitive than wild-type oocytes 
to exogenous DNA damage insults due to compromised HR 
repair. To understand the in vivo effect of OOEP depletion 
under the physiological conditions, we collected ovaries from 
4-week-old and 16-week-old wild-type and Ooep!"”"™ females 
and compared the number of oocytes within primordial and 
primary follicles, which express P63 and are subject to DNA 
damage-induced atresia (Suh et al., 2006). Intriguingly, at 
4-week old, the numbers of oocytes in primordial and primary 
follicles were significantly higher in Ooep*™™*™ females than in 
wild-type females. This probably reflected the developmental 
halt caused by the accumulated endogenous DNA damage in 
part of the oocytes. In contrast, at 16-week old, the numbers 
of oocytes in primordial and primary follicles were statistically 
lower in Ooep'"™ females than in wild-type females (Figure 
4B). Together, these results support that OOEP protects 
oocytes from endogenous or exogenous DNA DSB-induced 
apoptosis and preserves oocyte quantity under DNA damage 
insults. 

DNA damage in fully grown oocytes compromises the 
completion of meiosis by inducing a delay in meiosis resumption 
(germinal vesicle breakdown, GVBD) (Marangos & Carroll, 2012) 
and the activation of spindle assembly checkpoint (SAC) at 


metaphase | (Collins et al., 2015; Marangos et al., 2015). 


We investigated the influence of OOEP loss on the kinetics 
of GVBD and the first polar body extrusion (PBE) during 
oocyte in vitro maturation. Under normal conditions, the 
Ooep'"”™ oocytes showed a 4-h delay in reaching 50%-GVBD 
compared with wild-type oocytes. In addition, ~20% less 
oocytes completed GVBD in the absence of OOEP compared 
with the wild-type counterparts at the end of culture (Figure 
4C). Likewise, the Ooep""” oocytes displayed a 2-h delay in 
reaching 50%-PBE. At the end of culture, 20% less oocytes 
completed PBE compared with the wild-type counterparts 
(Figure 4D). Thus, OOEP was essential for efficient meiosis 
resumption and completion by preserving oocyte genomic 
integrity. Decline in DNA DSB repair competence is considered 
a major factor contributing to oocyte and ovarian aging at 
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advanced maternal age (Oktay et al., 2014, 2015; Titus et 
al., 2013; Xu et al., 2015). We examined if Ooep expression 
declined with maternal aging. mRNA expression levels of the 
Ooep gene measured in single GV oocytes of old females 
(40-week old) were lower than those from young females 
(8-week old) (Figure 4E). Thus, Ooep expression exhibited 
age-dependent decline and may contribute to ovarian failure 
and reproductive aging under physiological or pathological 
conditions. 


DISCUSSION 


DNA damage in oocytes can cause infertility, abortion, and 
birth defects, leading to reproductive failure. Among all 
types of DNA damage, DSBs are the most deleterious and 
can substantially impair genome integrity. Experimental and 
clinical studies on mice and humans have revealed that intact 
DNA DSB repair in oocytes is necessary to preserve oocyte 
quantity and quality. Decline in DNA DSB repair has a causal 
relationship with ovarian failure, menopause, and infertility 
(Oktay et al., 2014, 2015; Titus et al., 2013; Xu et al., 2015). 
Genome-wide association studies on menopause timing have 
also highlighted the strong association of early menopause 
with DNA damage response pathways where many genes 
are involved in DNA repair, particularly the HR pathway 
(e.g., BRCA1 pathway) and checkpoint response (cell cycle 
response and apoptosis) (Day et al., 2015; Perry et al., 2013). 
Thus, investigating how oocytes respond to DNA DSBs and 
what molecules are involved in HR-mediated DSB repair will 
shed light on oocyte and ovarian aging. In this study, we 
presented the following lines of evidence to support the notion 
that OOEP may participate in HR-mediated DNA DSB repair 
by regulating the ATM-RAD51 axis in oocytes. Firstly, after 
etoposide treatment, y-H2AX foci, which monitor DNA DSBs, 
were gradually resolved in wild-type oocytes but persisted 
in Ooep null oocytes, although this defect was rescued by 
re-expression of OOEP. Secondly, wild-type oocytes evoked 
the activation of ATM, a central coordinator in DNA DSB repair 
processes, in response to exogenous DNA damage. However, 
Ooep-/- oocytes failed to efficiently ignite ATM activation and 
this defect could be rescued by the re-introduction of OOEP 
proteins into mutant oocytes. Finally, RAD51, a recombinase 
essential for HR-mediated DNA DSB repair, was recruited to 
damage sites and formed foci in wild-type oocytes in response 
to etoposide treatment. In contrast, fewer RAD51 proteins were 
recruited to and formed foci at DNA DSBs in Ooep-/- oocytes, 
although this defect could also be rescued by exogenous 
OOEP. 

Due to the compromised DNA damage repair, Ooep-/- 
oocytes displayed an accumulation of endogenous DNA 
damage, as determined by the comet assay, an unambiguous 
method that measures the extent of DNA damage on a 
single cell basis. As DNA damage can induce atresia 
of oocytes in primordial and primary follicles through the 
p63-mediated pathway (Gonfloni et al., 2009; Suh et al., 
2006), we proposed that the accumulation of endogenous 
DNA damage in Ooep-/- oocytes may drive the atresia 


of oocytes and decrease the oocyte number in primordial 
and primary follicles under physiological conditions. Indeed, 
in young adult mice (16-week old), a mild but statistically 
significant decrease in oocyte quantity within the primordial 
and primary follicles was detected in the Ooep-/- ovaries 
compared to the wild-type counterparts. Oocytes in primordial 
and primary follicles from Ooep-/- females also displayed 
higher sensitivity to exogenous DNA insults and were prone 
to apoptosis. These observations suggest that OOEP is 
necessary for the preservation of oocyte quantity in normal and 
DNA-damaged conditions. OOEP depletion not only impairs 
the survival of oocytes within primordial and primary follicles, 
but also affects meiotic maturation and early embryonic 
development. Compared with wild-type oocytes, Ooep-/- GV 
oocytes exhibited delays in meiosis resumption as well as 
progression to the metaphase ll stage. Moreover, embryos 
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derived from Ooep-/- oocytes were arrested at the 2-cell stage 
(Li et al., 2008). This developmental arrest could be due 
to the DNA damage-induced cell cycle checkpoint (Xu et 
al., 2015). Thus, OOEP is essential for preserving oocyte 
quality. Interestingly, the mRNA expression of Ooep in mouse 
oocytes was reduced with advanced maternal age. Several 
key HR pathway genes, including Brca1, Mre11, Atm, and 
Rad51, also display age-dependent expression decrease in 
both human and mouse oocytes (Day et al., 2015; Oktay 
et al., 2015; Titus et al., 2013). Therefore, the coordinated 
decrease in the expression of Ooep and other canonical HR 
genes contributes to ovarian failure and reproductive aging 
under physiological and pathological conditions. OOEP is 
conservatively expressed in human oocytes (Zhu et al., 2015), 
implying that it may play similar roles in repairing DNA DSBs 
via the HR pathway in human oocytes. 
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Figure 4 Ooep protects oocytes from DNA DSB-induced apoptosis and meiosis delay 

A: Postnatal day 5 (P5) ovaries from WT and Ooep*”“*™ females were subjected to the same cisplatin treatment. TUNEL examination showed that the numbers 
of apoptotic germ cells (DDX4* TUNEL’, arrowheads) were higher in Ooep'”'” ovaries than in WT ovaries. B: HE staining revealed that ovaries from 4-week-old 
Ooep'™'” females contained more primordial and primary follicles than those from WT females (left panel). At 16-week old, Ooep!”'” ovaries contained less 
primordial and primary follicles than WT ovaries (right panel). Ten ovarian sections were examined for each animal and five females were used at each stage. C: 
The germinal vesicle breakdown (GVBD) of Ooep""”” oocytes was delayed compared to the WT oocytes under normal conditions. There was a 4-h difference 
in the timing of 50%-GVBD between the two samples. More than 80 oocytes were used in each sample. D: Kinetics of the first polar body extrusion (PBE) in 
WT and Ooep!”'” GV oocytes under normal conditions. Ooep""” oocytes displayed 2-h delay in reaching 50%-PBE. More than 50 oocytes were used in each 
sample. E: mRNA expression of Ooep declined with maternal aging in GV oocytes. Fourteen single GV oocytes from two females were examined at each stage. 


Data are presented as mean+ SEM. *: P<0.05, **: P<0.01, ***: P<0.001. 


OOEP proteins are distributed in cytoplasm and localized at 
the subcortex of mouse and human oocytes (Li et al., 2008; 
Zhu et al., 2015). OOEP was also excluded from the DNA 
DSB sites in oocytes treated with etoposide (data not shown). 
Thus, the functions of OOEP in regulating HR-mediated DNA 


DSB repair are likely indirect. OOEP is predicted to contain 
an atypical hnRNP K homology (KH) domain implicated in 
RNA binding (Wang et al., 2014). RNA binding proteins have 
been recognized as crucial players in the HR repair pathway 
(Adamson et al., 2012). The involvement of RBPs in regulating 
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ovarian aging has also been reported (Day et al., 2015). For 
instance, RNA binding protein FMRP encoded by fragile X 
mental retardation (Fmr1?) gene plays an essential role in 
preserving ovarian function in mice and humans (Ascano et al., 
2012; Lu et al., 2012). In the future, elucidating the molecular 
mechanism of OOEP is warranted to better understand the 
function of RBPs in preserving the genome integrity of germ 
cells and early embryos. 


CONCLUSIONS 


To the best of our knowledge, this work is the first to 
provide evidence that maternal OOEP may participate in 
HR-mediated DNA DSB repair by regulating the ATM-RAD51 
axis in oocytes. Thus, OOEP may participate in the regulation 
of genome stability in oocytes and contribute to ovarian failure 
and reproductive aging under physiological or pathological 
conditions. 
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ABSTRACT 


A new freshwater goby, Rhinogobius immaculatus sp. 


nov., is described here from the Qiantang River in 
China. It is distinguished from all congeners by the 
following combination of characters: second dorsal-fin 
rays |, 7-9; anal-fin rays I, 6-8; pectoral-fin rays 
14—15; longitudinal scales 29-31; transverse scales 
7-9; predorsal scales 2-5; vertebrae 27 (rarely 28); 
preopercular canal absent or with two pores; a red 
oblique stripe below eye in males; branchiostegal 
membrane mostly reddish-orange, with 3-6 irregular 
discrete or connected red blotches on posterior 
branchiostegal membrane and lower operculum in 
males; caudal-fin base with a median black spot; and 
no black blotch on anterior part of first dorsal fin in 
males. 


Keywords: Gobiidae; Rhinogobius; New species; 
Qiantang River; China 


INTRODUCTION 


The freshwater goby genus Rhinogobius Gill, 1859, is currently 
comprised of 74 valid species (Huang et al., 2016; Suzuki et 
al., 2017; Takahashi & Okazaki, 2017) widely distributed in 
East Asia, including Russia (Bogutskaya et al., 2008), Japan 
(Akihito et al., 2002), Korea (Regan, 1908b), China (Chen & 
Shao, 1996; Wu & Zhong, 2008), Philippines (Herre, 1927), 
Vietnam (Chen & Kottelat, 2005), Laos (Chen & Kottelat, 2003; 
Kottelat, 2001), Cambodia (Rainboth, 1996), and Thailand 
(Chen et al., 1999a). Most species of Rhinogobius from the 
islands of Japan and Taiwan are amphidromous (Chen & Shao, 
1996; Lee & Chang, 1996; Sakai et al., 2000), whereas most 
species from eastern continental Asia and Hainan Island are 
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non-diadromous (landlocked) (Chen et al., 1999a, 2002; Chen 
& Kottelat, 2005; Chen & Miller, 2014; Huang & Chen, 2007; Li 
& Zhong, 2009). 

In total, 44 species of Rhinogobius have been recorded in 
China (Chen et al., 2008; Chen & Miller, 2014; Huang et al., 
2016; Huang & Chen, 2007; Li et al., 2007; Li & Zhong, 2007, 
2009; Wu & Zhong, 2008; Yang et al., 2008), eight of which 
have been reported from the Qiantang River basin originating 
in southeastern Anhui Province to eastern Zhejiang Province. 
These species include R. aporus (Zhong € Wu, 1998), R. davidi 
(Sauvage € de Thiersant, 1874), R. cliffordpopei (Nichols, 
1925), R. leavelli (Herre, 1935a), R. lentiginis (Wu & Zheng, 
1985), R. niger Huang, Chen € Shao, 2016, R. similis Gill, 
1859, and A. wuyiensis Li & Zhong, 2007 (Chen et al., 1990; Li, 
2011; Liu et al., 2011; Suzuki et al., 2016; Zheng, 1989; Zheng 
& Wu, 1985). Herein, we describe a new species from three 
tributaries of the Qiantang River, China. 


MATERIALS AND METHODS 


Specimens for morphological examination were initially 
preserved in 6% formalin for seven days, and then transferred 
into 70% ethanol for permanent storage. Methods for 
morphometric measurements and meristic counts followed 
Nakabo (2002), with exceptions as indicated: Standard length 
(SL), head length, snout length, predorsal length, and preanal 
length were measured from the tip of the upper lip; Head depth 
and width were taken at the posterior margin of the preopercle; 
Body depth and width were taken at the origin of the anal fin. 
Vertebrae were counted from radiographs using the Kodak DXS 
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4000 system, and 3D reconstructed CT scans were made 
with the NSI-x50 system. Notations of cephalic sensory-canal 
pores and sensory-papillae rows followed Akihito et al. (2002) 
and Suzuki et al. (2017). Examined specimens in this study 
were deposited in the Biological Museum, Fudan University, 
Shanghai (FDU) and Shanghai Ocean University, Shanghai 
(SOU/SFU/SFC). 


RESULTS 


Rhinogobius immaculatus sp. nov. 
Figures 1-5; Table 1. 


Holotype: FDU 1010001, male, 25.7 mm SL; a tributary of 
Qiantang River, Fuchunjiang Town, Tonglu County, Zhejiang 
Province, China; 7 October 2010. 


Paratypes: FDU 1010002, female, 26.3 mm SL; same data as 
holotype. FDU 0905001—0905002, 2 males, 19.2-19.7 mm SL; 
FDU 0905003-0905006, 4 females, 20.9-21.4 mm SL; same 
locality as holotype; 4 May 2009. FDU 1107001-1107003, 3 
males, 20.5-22.1 mm SL; FDU 1107004—1107014, 11 females, 
22.2-25.2 mm SL; a tributary of Qiantang River, Xikou Town, 


Xiuning County, Anhui Province, China; 23 July 2011. FDU 
1710001-1710002, 2 females, 20.4—20.7 mm SL; a tributary of 
Qiantang River, Dongyangjiang Town, Dongyang City, Zhejiang 
Province, China; 3 October 2017. 


Diagnosis: Most similar to Rhinogobius wuyanlingensis in 
number of vertebrae (27) and preopercular canal pores (2 or 
O vs. 2), but differing by fewer pectoral-fin rays (14-15 vs. 
17-18), fewer anal-fin rays (l, 6-7 vs. |, 8), fewer transverse 
scales (7—9, modally 8 vs. 9-10), absence of a black blotch 
on anterior part of first dorsal fin in males (vs. present), 
and branchiostegal membrane mostly reddish-orange, with 
irregular blotches posteriorly in males (vs. with red stripes). 


Description: The morphometric and meristic data of the 
holotype and paratypes are shown in Table 1. The following 
features can describe the new species: Body cylindrical 
anteriorly and compressed posteriorly. Head sub-cylindrical. 
Eye large, dorsolateral. Mouth oblique, with lower jaw tip 
anterior-most, jaw forming 45° angel with body axis; corner of 
mouth reaching below anterior margin of orbit in adult males, 
not reaching below anterior margin in females. Vertebral counts 
11+16=27 (10) or 11+17=28 (1) (Figure 2). 





Figure 1 Rhinogobius immaculatus sp. nov. 


A: FDU 1010001, holotype, 25.7 mm SL, male; B: FDU 1107004, paratype, 24.4 mm SL, female. 
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Table 1 Meristic and morphometric data for holotype and paratypes of Rhinogobius immaculatus sp. nov. 





Holotype Paratypes 

Male Males Females 
Number 1 5 18 
Standard length (mm) 25:7 19.2-22.1 20.9-25.2 
First dorsal-fin rays VI VI (5) V (1); VI (15) 
Second dorsal-fin rays 1,9 1,7 (2); 1,8 (3) 1,7 (1); 1,8 (15); 1,9 (2) 
Anal-fin rays 1,7 1,6 (2); 1,7 (3) 1,6 (2); 1,7 (14); 1,8 (2) 
Pectoral-fin rays 15 15 (5) 14 (7); 15 (11) 
Longitudinal scales 30 29 (2); 30 (2); 31 (1) 29 (8); 30 (8); 31 (2) 
Transverse scales 7 8 (5) 7 (6); 8 (11); 9(1) 
Predorsal scales 4 2 (1); 3 (2);4(1);5(1) 2 (4); 3 (6); 4 (4); 5 (4) 
Vertebrae 27 unknown 27 (9); 28 (1); unknown (8) 
Morphometry 
% Standard length 
Head length 27.8 28.1-29.7 25.0-27.7 
Head depth 15.0 15.2-16.1 13.8-14.8 
Head width 19.4 17.9-20.6 17.1-18.3 
Body depth of anal-fin origin 14.1 14.2-15.1 14.1-16.3 
Body width of anal-fin origin 11.2 12.3-13.3 10.7-12.5 
Snout length 6.7 6.1-7.3 5.6-6.9 
Lower jaw length 8.6 7.2-9.3 6.3-7.6 
Orbit diameter 7.9 6.7—8.3 6.7—8.5 
Predorsal length 36.7 35.5-38.1 35.0-36.3 
Preanal length 61.1 58.4-59.7 58.7-61.9 
Caudal peduncle length 25.0 24.8-27.1 24.7-28.6 
Caudal peduncle depth 7.6 7.7-10.9 7.4-11.1 
Depressed 1st dorsal-fin length 19.2 17.5-18.3 16.8-18.8 
Depressed 2nd dorsal-fin length 32.7 27.0-31.0 26.5-29.7 
Depressed anal-fin length 24.0 22.4-25.5 21.6-24.2 
Pectoral-fin length 24.4 22.9-25.8 20.5-23.8 
Pelvic-fin length 16.7 17.0-18.1 15.3-18.7 
Caudal-fin length 24.7 23.0-26.5 21.6-25.3 


Numbers in parentheses are numbers of specimens with a given count. 





Figure 2 Rhinogobius immaculatus sp. nov. 
FDU 1107006, paratype, 22.2 mm SL, female. 3D reconstructed CT scan (with concealment of left skull). 
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First dorsal-fin rays V—VI (modally VI); second dorsal-fin 
rays |, 7-9 (modally |, 8); anal-fin rays l, 6-8 (modally |, 
7); pectoral-fin rays 14-15 (modally 15); pelvic-fin rays I, 5; 
segmented caudal-fin rays 9+8, including branched rays 7+7; 
dorsal procurrent rays 6-8, ventral procurrent rays 5—7. First 
dorsal fin with third or fourth spine longest, no filamentous 
spines; rear tip not reaching origin of second dorsal fin 
when depressed in both sexes. Second dorsal and anal 
fins short-based, tip of depressed rays far from dorsal and 
ventral origins of procurrent caudal-fin rays. Origin of anal fin 
inserted below base of third and fourth rays of second dorsal 
fin. Pectoral fin elliptical, central rays longest; rear extension far 


from vertical of anus when depressed. Pelvic fin disc rounded. 


Rear edge of caudal fin rounded. 
Longitudinal scales 29-31 (28-30 on body, 0—1 on caudal 


fin); transverse scales 7—9 (modally 8); predorsal scales 2-5. 


Body with moderately ctenoid scales. Anterior predorsal area, 
head, pectoral base, and prepelvic area naked. Posterior 
predorsal area and belly with cycloid scales. Anterior-most 
predorsal scale not reaching vertical through upper end of gill 
opening. 

Head pores present. Nasal extension of anterior 
oculoscapular canal with terminal pores B’ at vertical between 
anterior and posterior nostrils. Anterior interorbital section 
of oculoscapular canal separated, with paired pore C. Single 
pore D in posterior interorbital region. Postorbital region 
with paired pore E. Lateral section of anterior oculoscapular 
canal with anterior pore F and terminal pore H’. Posterior 
oculoscapular canal short, reduced (absent in 14 specimens), 
with two terminal pores K’ and L. Gap between anterior and 
posterior oculoscapular canals larger than length of posterior 
oculoscapular canal. Preopercular canal short, reduced 
(absent in eight specimens), with two terminal pores M’ and 
O’. Sensory-papillae row a short, with two or three papillae 
below orbit. Row b short, about half length of orbit. Rows c 
and d longer, not extending to vertical line of rear margin of 
orbit. Single cp papilla. Row f paired (Figure 3). 


Color in life: Ground color light brown. Snout with pair of 
reddish brown stripes united at tip of snout. A reddish oblique 
stripe below eye in males, not reaching rear edge of mouth; 
obscure in females. Cheek and opercle with irregular reddish 
lines, branchiostegal membrane mostly reddish-orange, with 
3—6 irregular discrete or connected red blotches on posterior 
branchiostegal membrane and lower operculum in males; 
absent in females. Flank with 6-7 irregular discrete or 
connected black blotches. First dorsal fin with 3-4 rows of 
interphased black and white spots. A black blotch on anterior 
part of first dorsal fin absent in both sexes. Second dorsal fin 
with 4—5 rows of interphased black and white spots. Pectoral 
fin proximally white, posterior part with 3—4 rows of interphased 
black and white spots. Pectoral-fin base with irregular blackish 
pigmentation, usually darker in upper part. Caudal fin with 5-6 
rows of interphased black and white spots. Caudal-fin base 
with a median black blotch. Pelvic fin and anal fin with slight 
irregular black pigmentation (Figures 4A-—B, 5). 


Distribution and ecology: Known only from streams of the 
Qiantang River basin in Zhejiang and Anhui Provinces, China 
(Figure 6). Most often found in shallow (10-50 cm deep) 
low-gradient streams, with sand and gravel mixed substrate. 

Adult Rhinogobius immaculatus sp. nov. are small in size. 
The smallest female with mature oocytes was 22.4 mm SL. The 
largest specimen collected in the field was 26.3 mm SL. The 
largest captive specimen kept in an aquarium for 29 months 
was 32.8 mm SL. 


Etymology: The specific name, immaculatus, is derived 
from Latin in (without) and maculatus (spotted), an adjective, 
alluding to the absence of a black blotch on the anterior part of 
the first dorsal fin in adult males. 


DISCUSSION 


Number of vertebrae is frequently used for species 
identification in the genus Rhinogobius (Chen et al., 2008; 
Huang et al., 2016; Lee & Chang, 1996; Suzuki et al., 2017; 
Takahashi & Okazaki, 2017; Yang et al., 2008). Among the 
current 74 valid species in Rhinogobius, 36 possess 27 or 
more vertebrae (Huang et al., 2016), as also found in the newly 
described species, and 28 possess less than 27 vertebrae 
(Suzuki et al., 2016, 2017; Takahashi & Okazaki, 2017). The 
vertebral number of the remaining species remains unknown. 
In the genus Rhinogobius, only R. wuyanlingensis Yang, Wu 
8 Chen, 2008 and R. lindbergi Berg, 1933 are known to have 
two preopercular canal pores (Huang et al., 2016; Sakai et al., 
2000). Rhinogobius immaculatus sp. nov. closely resembles 
R. wuyanlingensis in number of vertebrae (27), preopercular 
canal pores (2 or O vs. 2), presence of predorsal scales, and 
small adult size, but differs from R. wuyanlingensis in fewer 
pectoral-fin rays (14-15 vs. 17-18), fewer anal-fin rays (I, 
6-7 vs. |, 8), fewer transverse scales (7—9, modally 8 vs. 
9-10), absence of a black blotch on the anterior part of the first 
dorsal fin in males (vs. present), and branchiostegal membrane 
mostly reddish-orange, with irregular blotches posteriorly in 
males (vs. with red stripes) (Yang et al., 2008). Rhinogobius 
immaculatus sp. nov. shares the same number of vertebrae 
(27) and preopercular canal pores (2 or O vs. 2) with 
R. lindbergi, but differs in fewer pectoral-fin rays (14-15 vs. 
19-20), fewer anal-fin rays (l, 6-7 vs. |, 9), and more predorsal 
scales (2-5 vs. 0) (Berg, 1933; Sakai et al., 2000). 
Rhinogobius immaculatus sp. nov. can be distinguished 
from the other 34 species with 27 or more vertebrae as 
follows: from R. albimaculatus Chen, Kottelat & Miller, 1999a, 
R. boa Chen & Kottelat, 2005, R. changtinensis Huang & Chen, 
2007, R. chiengmaiensis Fowler, 1934, R. duospilus (Herre, 
1935b), R. filamentosus (Wu, 1939), R. flumineus (Mizuno, 
1960), R. henryi (Herre, 1938), R. honghensis Chen, Yang & 
Chen, 1999c, R. lineatus Chen, Kottelat 8 Miller, 1999a, R. 
linshuiensis Chen, Miller, Wu & Fang, 2002, R. longyanensis 
Chen, Cheng & Shao, 2008, R. lungwoensis Huang € Chen, 
2007, R. maculicervix Chen € Kottelat, 2000, R. mekongianus 
(Pellegrin & Fang, 1940), R. milleri Chen € Kottelat, 2003, R. 
nammaensis Chen & Kottelat, 2003, R. niger, R. ponkouensis 
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Huang & Chen, 2007, R. sulcatus Chen € Kottelat, 2005, R. 
taenigena Chen, Kottelat & Miller, 1999a, R. vermiculatus Chen 
& Kottelat, 2003, R. wangchuangensis Chen, Miller, Wu & Fang, 
2002, R. wangi Chen & Fang, 2006, R. xianshuiensis Chen, 
Wu & Shao, 1999b, and R. yaoshanensis (Luo, 1989) by fewer 
preopercular canal pores (2 or 0 vs. 3) and absence of a black 


blotch on the anterior part of the first dorsal fin in males (vs. 


present) (Chen et al., 2008; Huang et al., 2016; Wu & Zhong, 
2008); from R. genanematus Zhong € Tzeng, 1998, and R. 
parvus (Luo, 1989) by more predorsal scales (2-5 vs. 0-4, 
usually 0 in R. parvus, and 0 in R. genanematus) and fewer 
preopercular canal pores (2 or 0 vs. 3) (Chen et al., 2008); from 


R. cheni (Nichols, 1931) by fewer longitudinal scales (29-31 vs. 


34) and fewer preopercular canal pores (2 or 0 vs. 3) (Chen et 
al., 2008); from R. szechuanensis (Tchang, 1939) by presence 
of oculoscapular canal (vs. absent) and more predorsal scales 
(2-5 vs. 0) (Chen et al., 2008; Wu 8 Zhong, 2008); from R. 
davidi and R. multimaculatus (Wu 8 Zheng, 1985) by more 
predorsal scales (2-5 vs. 0) and fewer vertebrae (27 vs. 28 in 
R. davidi, 29 in R. multimaculatus) (Chen & Miller, 1998; Yang 
et al., 2008; Zheng & Wu, 1985); from R. rubromaculatus Lee 
& Chang, 1996 by fewer transverse scales (7-9 vs. 10-13) 
and fewer predorsal scales (2-5 vs. 9-13) (Chen € Shao, 


1996); and from R. lentiginis by fewer transverse scales (7—9 vs. 


10-11), absence of spots on cheek (vs. present), and absence 
of a black blotch on the anterior part of the first dorsal fin in 
males (vs. present) (Li, 2011; Zheng € Wu, 1985). 

In addition to differences in vertebral number, Rhinogobius 
immaculatus sp. nov. can be distinguished from the 28 
species with less than 27 vertebrae as follows: from R. 
aporus, R. changjiangensis Chen, Miller, Wu & Fang, 2002, 
R. leavelli, R. nandujiangensis Chen, Miller, Wu & Fang, 2002, 
R. reticulatus Li, Zhong & Wu, 2007, R. rubrolineatus Chen & 
Miller, 2008, R. sagittus Chen € Miller, 2008, R. sangenloensis 
Chen & Miller, 2014, R. variolatus Chen & Kottelat, 2005, R. 
virgigena Chen € Kottelat, 2005, and R. wuyiensis by absence 
of a black blotch on the anterior part of the first dorsal fin in 
males (vs. present) (Herre, 1935a; Li & Zhong, 2007; Wu & 
Zhong, 2008; Zhong & Wu, 1998); from R. biwaensis Takahashi 
& Okazaki, 2017, R. brunneus (Temminck & Schlegel, 1845), 
R. candidianus (Regan, 1908a), R. delicatus Chen & Shao, 
1996, R. fluviatilis Tanaka, 1925, R. formosanus Oshima, 1919, 
R. gigas Aonuma & Chen, 1996, R. henchuenensis Chen & 
Shao, 1996, R. kurodai (Tanaka, 1908), R. lanyuensis Chen, 
Miller & Fang, 1998, R. maculafasciatus Chen € Shao, 1996, 
R. mizunoi Suzuki, Shibukawa € Aizawa, 2017, R. nagoyae 
Jordan & Seale, 1906, R. nantaiensis Aonuma & Chen, 1996, 
and R. ogasawaraensis Suzuki, Chen & Senou, 2012 by fewer 
pectoral-fin rays (14-15 vs. more than 16) and presence of a 
median black blotch on the caudal fin base (vs. absent) (Akihito 
et al., 2002; Chen & Shao, 1996; Suzuki & Chen, 2011; van 
Oijen et al., 2011); from R. similis by fewer pectoral-fin rays 
(14-15 vs. 18-19) and fewer predorsal scales (2-5 vs. 8-12) 
(Gill, 1859; Suzuki et al., 2016); and from R. zhoui Li & Zhong, 
2009 by fewer predorsal scales (2-5 vs. 10-12) and absence 
of white margin on each median fin (vs. present). 
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Rhinogobius immaculatus sp. nov. can be distinguished 
from the 10 species with unknown vertebral number as follows: 
from R. bedfordi (Regan, 1908b), R. bucculentus (Herre, 1927), 
R. carpenteri Seale, 1910, and R. philippinus (Herre, 1927) by 
fewer longitudinal scales (29-31 vs. more than 35); from R. 
cliffordoopei by more predorsal scales (2-5 vs. 0) (Nichols, 
1925); from R. fukushimai Mori, 1934, and R. imfasciocaudatus 
Nguyen & Vo, 2005 by fewer pectoral-fin rays (14-15 vs. 18) 
(Nguyen, 2005); from R. liui Chen € Wu, 2008 by fewer 
pectoral-fin rays (14-15 vs. 19), fewer longitudinal scales 
(29-31 vs. 36-39), and more predorsal scales (2-5 vs. 0) (Wu 
& Zhong, 2008); from R. shennongensis (Yang € Xie, 1983) by 
fewer pectoral-fin rays (14-15 vs. 18-19) and fewer longitudinal 
scales (29-31 vs. 31-33); and from R. sowerbyi Ginsburg, 
1917 by fewer longitudinal scales (29-31 vs. 35-36) and fewer 
anal-fin rays (l, 6-7 vs. 1,8). 


COMPARATIVE MATERIAL 


Rhinogobius aporus: FDU 200910103, topotypes, 14 
specimens, 23.0-31.5 mm SL; Ou River, China: Zhejiang 
Province: Jinyun County; 9 October 2009. FDU 201010106, 16 
specimens, 21.9-37.5 mm SL; Qiantang River, China: Zhejiang 
Province: Wuyi County; 9 October 2010. FDU 201710101, 
23 specimens, 23.5-35.7 mm SL; Ou River, China: Zhejiang 
Province: Pan’an County; 4 October 2017. 

Rhinogobius cliffordpopei: FDU 201010105, 6 specimens, 
25.3-32.0 mm SL; Qiantang River, China: Zhejiang Province: 
Yongkang City; 8 October 2010. FDU 201309101, 11 
specimens, 24.2-36.9 mm SL; Yangtze River, China: Hunan 
Province: Lianyuan City; 26 September 2013. 

Rhinogobius davidi: FDU 200905101, 6 specimens, 
22.5-27.0 mm SL; Qiantang River, China: Zhejiang Province: 
Tonglu County; 24 May 2009. FDU 201110101, 18 specimens, 
18.4-37.5 mm SL; Qiantang River, China: Zhejiang Province: 
Kaihua County; 7 October 2011. 

Rhinogobius changtinensis: FDU 201711101, topotypes, 
18 specimens, 28.2-38.9 mm SL; Han River, China: Fujian 
Province: Changting County; 28 November 2017. 

Rhinogobius duospilus: FDU 200711101, 24 specimens, 
21.8-33.0 mm SL; China: Guangdong Province: Shenzhen 
City; November 2007. FDU 201310101, 6 specimens, 
31.1-39.6 mm SL; Pearl River, China: Guangxi Province: 
Jinxiu County; 1 October 2013. FDU 201712101, 6 specimens, 
22.2-28.7 mm SL; Pearl River, China: Guangxi Province: 
Nanning City; 26 December 2017. 

Rhinogobius formosanus: FDU 201512101, 13 specimens, 
29.8-54.2 mm SL; China: Fujian Province: Fuging City; 12 
December 2015. 

Rhinogobius genanematus: FDU 200910102, topotypes, 11 
specimens, 19.5-28.4 mm SL; Ling River, China: Zhejiang 
Province: Xianju County; 8 October 2009. FDU 201007102, 
14 specimens, 22.6-32.3 mm SL; Ling River, China: Zhejiang 
Province: Tiantai County; 27 July 2010. 

Rhinogobius gigas: FDU 201409101, 4 specimens, 
40.8-62.4 mm SL; Hsiukuluan River, China: Taiwan: Taitung 
County; 20 September 2014. 


Rhinogobius honghensis: FDU 201002101, 28 specimens, 
35.8-49.4 mm SL; Red River, China: Guangxi Province: Napo 
County; February 2010. 

Rhinogobius leavelli: FDU 201010104, 7 specimens, 
34.5-44.2 mm SL; Qiantang River, China: Zhejiang Province: 
Tonglu County; 7 October 2010. FDU 201108101, 9 
specimens, 23.1-51.3 mm SL; Lingshui River, China: Hainan 
Province: Baoting County; 21 August 2011. 

Rhinogobius lentiginis: FDU 201007102, topotypes, 17 
specimens, 25.4-40.3 mm SL; Ling River, China: Zhejiang 
Province: Tiantai County; 27 July 2010. FDU 200706101, 
21 specimens, 24.8-37.2 mm SL; Cao'e River (tributary of 
Qiantang River), China: Zhejiang Province: Xinchang City; 10 
June 2007. 

Rhinogobius lindbergi: FDU 201208101, topotypes, 5 
specimens, 19.5-25.1 mm SL; China: Heilongjiang Province: 
Harbin City; 28 August 2012. 

Rhinogobius linshuiensis: FDU 201108101, topotypes, 13 
specimens, 20.2-37.5 mm SL; Lingshui River, China: Hainan 
Province: Baoting County; 5 August 2011. 

Rhinogobius multimaculatus: FDU 201107101, topotypes, 
5 specimen, 19.7-35.2 mm SL; Tiaoxi River China: Zhejiang 
Province: Anji County; 24 July 2011. 

Rhinogobius niger: FDU 200910102, topotypes, 3 
specimens, 29.9-34.2 mm SL; Ling River, China: Zhejiang 
Province: Xianju County; 8 October 2009. FDU 201007101, 
15 specimens, 22.1-43.9 mm SL; Yong River, China: Zhejiang 
Province: Fenghua City; 26 July 2010. FDU 201007102, 10 
specimens, 17.5-33.7 mm SL; Cao’e River (tributary of Qiantang 
River), China: Zhejiang Province: Xinchang City; 27 July 2010. 

Rhinogobius reticulatus: SFU 07001, holotype, male, 33.6 
mm SL; Min River, China: Fujian Province: Fuzhou City; 
August 2005. SFU 07002, paratype, female, 30.4 mm SL; SFU 
07003-07006, paratypes, 4 males, 27.3-35.8 mm SL; SFU 
07007-07011, paratypes, 6 females, 25.6-34.5 mm SL; same 
data as holotype. FDU 201010102, 14 specimens, 20.5-27.8 
mm SL; Min River, China: Fujian Province: Pucheng County; 5 
October 2010. 

Rhinogobius rubromaculatus: FDU 201503101, 8 
specimens, 19.5-34.8 mm SL; Shihwen River, China: Taiwan: 
Pingtung County; 29 March 2015. 

Rhinogobius shennongensis: FDU 201404101, topotypes, 
5 specimens, 46.5-50.3 mm SL; Yangtze River, China: Hubei 
Province: Shennongjia Forestry District; 23 April 2014. 

Rhinogobius similis: FDU 201009101, 6 specimens, 
44.1-62.5 mm SL; Yangtze River, China: Jiangxi Province: 
Wuyuan County; 3 Sep. 2010. FDU 201010101, 5 specimens, 
39.5-52.4 mm SL; Qiantang River, China: Zhejiang Province: 
Changshan; County; 3 September 2010. FDU 201708101, 2 
specimens, 36.5-36.9 mm SL; Perfume River, Vietnam: Hue 
City; 29 August 2017. 

Rhinogobius szechuanensis: FDU 201005101, topotypes, 7 
specimens, 41.6-58.9 mm SL; Yangtze River, China: Sichuan 
Province: Qionglai City; May 2010. 

Rhinogobius yaoshanensis: FDU 201310101, topotypes, 6 
specimens, 23.6-41.5 mm SL; Pearl River, China: Guangxi 
Province: Jinxiu County; 1 October 2013. FDU 201310102, 


3 specimens, 39.6-43.9 mm SL; Pearl River, China: Guangxi 
Province: Bama County; 7 October 2013. 

Rhinogobius  wuyanlingensis: FDU 200710101, 8 
specimens, 176-27.7 mm SL; Ao River, China: Zhejiang 
Province: Pingyang City; 12 October 2007. FDU 200910101, 
19 specimens, 20.3-27.8 mm SL; Dajing River, China: 
Zhejiang Province: Yueging City; 4 October 2009. 

Rhinogobius wuyiensis: SFU 07101, holotype, male, 39.2 
mm SL; Qiantang River, China: Zhejiang Province: Wuyi 
County; July 2006. SFU 07102, paratype, female, 37.0 mm 
SL; SFU 07103-07109, paratypes, 7 males, 32.0-41.6 mm 
SL; SFU 07110-07118, paratypes, 8 females, 31.6-38.6 mm 
SL; same data as holotype. FDU 201010103, 8 specimens, 
19.9-27.5 mm SL; Qiantang River, China: Zhejiang Province: 
Lanxi City; 6 October 2010. FDU 201007101, 12 specimens, 
17.6-35.9 mm SL; Yong River, China: Zhejiang Province: 
Fenghua City; 26 July 2010. 





Figure 3 Cephalic lateral-line system of Rhinogobius 
immaculatus sp. nov. 
FDU 1107003, paratype, 20.5 mm SL. 
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Figure 4 Male (A) and female (B) Rhinogobius immaculatus sp. nov. 


Specimens collected from Xiuning County, Anhui Province, China. Specimens not preserved. 





Figure 5 Coloration of head and first dorsal fin of male Rhinogobius immaculatus sp. nov. 


Specimen collected from Xiuning County, Anhui Province, China. Specimen not preserved. 
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Figure 6 Sampling localities of Rhinogobius immaculatus sp. nov. 


x: type locality; a: other sampling localities. 


Rhinogobius xianshuiensis: SFC 2257-2258, paratypes, 2 
specimens, 24.8-28.8 mm SL; Xianshui River, China: Fujian 
Province: Xianyou County; 19 August 1994. FDU 201407101, 
topotypes, 9 specimens, 23.1-32.0 mm SL; Xianshui River, 
China: Fujian Province: Xianyou County; 16 July 2014. 

Rhinogobius zhoui: SOU 0804001, holotype, male, 33.4 
mm SL; Huang River, China: Guangdong Province: Lianhua 
County; April 2008. SOU 0804002, paratype, female, 31.7 mm 
SL; SOU 0804003-0804008, paratypes, 6 males, 26.6-36.1 
mm SL; SOU 0804009-0804013, paratypes, 5 females, 
30.3-32.1 mm SL; same data as holotype. 
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ABSTRACT 


Understanding the spatial patterns of human-wildlife 
conflict is essential to inform management decisions 
to encourage coexistence, but it is constrained by 
the lack of spatially-explicit data. We collected 
spatially-implicit data of human-wildlife conflicts from 
2009-2015 around Daxueshan Nature Reserve, 
Yunnan, China, and investigated the patterns and 
drivers of these conflicts. A questionnaire was also 
designed to capture local resident attitudes toward 
insurance-based compensation for the losses caused 
by targeted wildlife. We found that the Asiatic black 
bear (Ursus thibetanus) was the most conflict-prone 
animal around the reserve, followed by the rhesus 
macaque (Macaca mulatta) and Southeast Asian 
sambar (Cervus equinus). Conflicts were unevenly 
distributed among seasons, villages, and communities, 


with several grids identified as conflict hotspots. 


Poisson models revealed that human-bear conflicts 
were negatively related to distance to the reserve 
and proportion of forest, but positively correlated 
to the proportion of cropland. Binomial models 
showed that communities affected by crop depredation 
were positively correlated with the proportion of 
cropland and negatively correlated with distance 
to the reserve, whereas communities affected by 
livestock depredation were negatively correlated with 
the proportion of cropland. The insurance-based 
scheme has compensated over 90% of losses, to 
the satisfaction of 90.6% of respondents. Our results 
suggest that human-bear conflict could be potentially 
reduced by eliminating food crops near the reserve 


boundary and livestock grazing at conflict hotspots. 
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In addition, the insurance-based scheme could be 
replicated at a broader scale with improvement in loss 
assessment. 


Keywords: Human-wildlife conflict; Asiatic black bear; 
Spatial heterogeneity; Insurance scheme; Daxueshan 
Nature Reserve 


INTRODUCTION 


Human-wildlife conflict, which is defined as any interactions 
leading to negative impacts on the humans or wildlife involved 
(Pettigrew et al., 2012), is a worldwide conservation issue. 
Large- and medium-sized mammals are the primary animals 
of concern and include species such as lions (Panthera leo), 
snow leopards (Panthera uncia), and Asian elephants (Elephas 
maximus) (Bagchi & Mishra, 2006; Chen et al., 2013; Distefano, 
2005; Maclennan et al., 2009), which are highly valued by 
international tourists and researchers. However, local residents 
often experience variously negative impact due to the presence 
of wildlife, including crop raiding, livestock depredation, and 
human casualties (Dickman et al., 2011). In turn, local resident 
hostility can increase, thereby threatening wildlife conservation 
(Dickman, 2010). 

Prevention before conflict and mitigation after conflict are 
two general strategies used to tackle human-wildlife conflict 
(Marchal & Hill, 2009; Mishra et al., 2003; Pettigrew et 
al., 2012). Prevention is widely recognized as the better 
strategy (Goodrich, 2010; Treves & Karanth, 2003), and 
includes guarding and fencing of livestock, zoning of land, 
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and increasing prey abundance for carnivores (Guo et al., 
2012; Mishra et al., 2003). Determining the areas in which 
to implement such measures is the first concern, and thus 
understanding the spatial heterogeneity of conflict is the 
cornerstone to achieve cost-efficient prevention (Dickman et al., 
2011). However, studies on spatial heterogeneity are limited 
by the lack of long-term spatially-explicit conflict data, which 
require extensive labor to collect and standard approaches to 
assess (Pettigrew et al., 2012). Therefore, conflict mitigation 
measures are still widely used (Li et al., 2013; Mishra et al., 2003), 
among which monetary compensation is the most common. As 
a result, spatially-implicit conflict data are well documented by 
compensation schemes, which could allow promising insight into 


the spatial patterns of human-wildlife conflict (Chen et al., 2016b). 


In China, compensation schemes are funded by several 
levels of government, whereby monetary compensation is 
offered for the losses caused by targeted animals listed as 


nationally protected species under the National Wildlife Law. 


Within Yunnan, the first scheme was funded by the provincial 
government from 1999-2008 and the second scheme was 
funded by the state government from 2009-2013. To date, 
the insurance-based scheme is the most expansive, and 
was initiated in 2007 to mitigate human-elephant conflict at 
Xishuangbanna, Yunnan, China. In this scheme, adjustors 
from the insurance company assessed the losses and recorded 
the conflict data, including the household affected, date of 
conflict, wildlife involved, type of damage, conflict locality, loss 
assessment, and compensation amount. Previous studies 
have reported that a high proportion of local residents were 
dissatisfied with this scheme after suffering significant losses 
caused by Asian elephant in south Yunnan and large-sized 
carnivore in the Qinghai-Tibetan Plateau (Chen et al., 2013, 
Chen et al., 2016a). However, local resident attitudes toward 
this scheme with different conflict-prone animals as well as 
conflict severity remain poorly evaluated. 

Here, we aimed to (1) identify the hotspots of human-wildlife 
conflict using spatially-implicit compensation data, (2) 
investigate factors potentially affecting the patterns of 
human-wildlife conflict, and (3) evaluate the attitudes of local 
residents toward insurance-based compensation schemes. 


MATERIALS AND METHODS 


Study area 

This study was conducted around the Daxueshan Nature 
Reserve (DNR), Yunnan, China (175.41 km?, E99°32’—-99° 43’ 
and N24°00’—24°12’). The elevation within the reserve 
ranges from 960 m to 3 504 m a.s.l., with diverse vegetation 
ranging from tropical seasonal rainforest to temperate alpine 
forest. Large- and medium-sized mammals are found within 
DNR, including the western black-crested gibbon (Nomascus 
concolor), northern pig-tailed macaque (Macaca leonina), 
rhesus macaque (Macaca mulatta), Southeast Asian sambar 
(Cervus equinus), Indian muntjac (Muntiacus vaginalis), wild 
boar (Sus scrofa), and Asiatic black bear (Ursus thibetanus) 
(Guo, 2006). In the area surrounding DNR, the predominant 
income of residents comes from cash crops (e.g., walnut, 


tobacco, and sugarcane) and livestock (e.g., goat and cattle). 


Data collection 

We collected human-wildlife conflict compensation data recorded 
by the DNR and insurance company from January 2009 
to September 2015, which consisted of 1 637 compensated 
conflicts. We assigned conflicts during 2014-2015 into 
1-km? grids by interviewing households in communities with 
the highest losses. High-resolution landscape images were 
loaded into Google Earth, then overlaid with the 1-km? grids 
labeled with localities. We also designed a questionnaire 
to collect information on local resident perceptions toward 
the insurance-based scheme (Supplementary Table S1). 
We collected socioeconomic information during 2009-2014 
from the Digital Village of Yunnan (http://www.ynszxc.gov.cn/) 
established by the Yunnan Provincial Government. 

We visually interpreted 360 points on the high-resolution 
imagery to develop a land-use map (Wang et al., 2016). We 
classified land use into five categories: forest, shrubland, 
cropland, construction site, and water body. 


Land use classification 

Two Landsat 8 OLI_TIRS images captured in the dry season 
in 2015 were used to develop the land-use map. We 
pre-processed the Landsat images as per Young et al. (2017). 
In addition to the spectral information provided by image bands, 
ancillary data were added to improve classification accuracy, 
including the Normalized Difference Vegetation Index (NDVI), 
NDVI texture, Digital Elevation Model (DEM), and slope. In 
total, 240 training points were used to perform supervised 
classification using a random forest algorithm with the help 
of the RStoolbox package in R (Leutner et al., 2017; R 
Development Core Team, 2011). The overall classification 
accuracy was 0.94 validated by 120 points. 


Statistical analysis 

Pre-analysis showed that loss was highly correlated with 
the number of conflicts (Bear-Livestock: r?=0.89, P<0.01; 
Bear-Crop; r?=0.78, P<0.01; Macaque-Crop: r?=0.96, P<0.01; 
Sambar-Crop: r?=1, P<0.01); hence, both loss and number of 
conflicts could be used as indicators of human-wildlife conflict. 
We excluded eight human-injury events from further analysis as 
they were stochastic and disproportionally accounted for 23.9% 
of losses. The predation rate of livestock was defined as the 
number killed divided by fatstock each year. 

We performed Chi-square tests of independence to 
examine whether human-wildlife conflict varied among villages, 
communities, and seasons. We mapped spatial heterogeneity 
of the losses from 2009 to 2014 with log transformation at the 
village scale by spatial Kriging interpolation (Chen et al., 2013). 
The villages were roughly represented by the location of the 
village committees. We assigned data into 1-km? grids from 
2014 to 2015 to identify conflict hotspots at the grid scale. 

We postulated that five factors potentially affected the 
number of human-bear conflicts (divided into crop depredation 
and livestock depredation) at the community level from 2009 
to 2014, including distance to DNR, forest, shrubland, cropland, 
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and construction site (as proportions within a 2-km radius of 
each community). Kendall tests confirmed that correlation 
of the five factors was lower than 0.7. Due to the 
large variation in the positive count data, we fitted the 
data with the “hurdle” model, which combined logistic 
regression for binomial responses (affected or not affected) 
and zero-truncated Poisson model for count data (Zeileis 
et al., 2008). Univariate analysis was carried out for the 
explanatory covariates, followed by multivariate analysis for 
significant covariates (excluding shrubland and construction 
sites, P<0.05) (Songhurst & Coulson, 2014). Moran's | test 
showed no significant spatial autocorrelation in the baseline 
model. The resulting models were ranked using Akaike's 
information criterion corrected for small sample sizes (AlCc) 
with the MuMIn package in R (Burnham & Anderson, 2002). 


RESULTS 
Description of human-wildlife conflict 


The 1 526 conflicts caused total economic losses of 
USD206 341. Asiatic black bear, rhesus macaque, and Southeast 
Asian sambar were the major animals involved (Table 1), 
contributing to 98.7% of the conflicts. The Asiatic black bear 
caused the greatest damage, accounting for 77.7% of conflicts 
and 88.2% of losses, followed by the rhesus macaque and 
Southeast Asian sambar (17.0% of conflicts and 8.6% of 
losses and 4.7% of conflicts and 1.9% of losses, respectively). 
For conflicts attributed to the Asiatic black bear, livestock 
depredation caused disproportional losses compared to the 
frequency of conflicts (49.6% vs. 18.2%), with an average 
predation rate of 2.57%. 


Table 1 Summary of human-wildlife conflicts during 2009-2014 around Daxueshan Nature Reserve 


Animal Frequency Economic loss (US$) 


Number of conflicts (n) 





Maize Goat Hive Bean Radish Buckwheat Other 





Ursus thibetanus 1185 181 999 
Macaca mulatta 260 17 736 
Cervus equinus 71 3 941 
Psittacula finschii 6 586 
Raptors 4 2079 
Sum in total 1 526 206 341 


Spatiotemporal pattern of human-wildlife conflict 

Conflicts varied among villages and communities (df=15, 
y?=67.779, P<0.01; df=76, y?=231.41, P<0.01) and grids 
(Figures 1, 2). We assigned villages into Parts 1 to 4. Part 
1 was comprised of four villages and accounted for 60.5% of 
losses, with 99.9% caused by the Asiatic black bear, such that 
the predation rate (9.5%) was much higher than the average. 
Part 2 was comprised of one village and accounted for 6.5% of 
losses, with 92.5% attributed to the Asiatic black bear. Part 3 
and Part 4 accounted for 9.7% and 10.3% of losses (with 80.6% 
and 58.6% attributed to the Asiatic black bear), respectively. 
At the grid scale, 92.6% of conflicts (113 out of 122) were 
located and assigned to 18 grids (Figure 2). Grids near 
communities Baishudi (BSD), Taoshudi (TSD), and Luchang 
(LC) suffered considerably from goat depredation (90.1% of 
losses). The grids near community Shijiaoyan (SJY) were 
inhabited by a group of rhesus macaques, which contributed to 
99.1% of conflicts. Crop depredation by the Asiatic black bear 
accounted for the main losses at the grids near communities 
Kuzhupeng (KZP) (100%) and Xinzhai (XZ) and Huomaoshu 
(HMS) (95.2%). 

Of the three predictors (Table 2), the probability of the 
community suffering cropland depredation was significantly 
negatively correlated with distance to DNR and positively 
correlated with the proportion of cropland, whereas the 
probability of the community suffering livestock depredation 
was significantly negatively correlated to the proportion of 
cropland. In the zero-truncated model with positive count data, 
the number of conflicts, regardless of type of depredation, was 
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significantly negatively correlated with the distance to DNR and 
proportion of forest, but significantly positively correlated to the 
proportion of cropland. 


Table 2 Coefficients of predictors fitted to the probability of 
the community being affected and number of conflicts 





: Livestock Livestock Crop Crop 
Covariates ; : p : ‘ f 
(Binomial) (Poisson) (Binomial) (Poisson) 
Distance to DNR z —0.77*** —0.95* —0.31*** 
Proportion of cropland -—0.72* 0. 39*** 1.04** 0.14*** 
Proportion of forest pe —0.42*** és —0.20*** 


*: P<0.05; **: P<0.01; ***: P<0.001; ~: P>0.05. 


The number of conflicts varied among seasons (df=3, 
y?=19.65, P<0.01) (Figure 3, 4). The Asiatic black bear 
predated goats most frequently from July to October (83.5% 
of losses) and 99.1% of maize-crop raiding by the black bear 
occurred from July to November. Rhesus macaques mainly 
caused damage to maize (62.3% of conflicts) and bean crops 
(33.5% of conflicts) between July and October. The Southeast 
Asian sambar caused damage to radish and buckwheat crops, 
with most of conflicts (78.9%) occurring from October to 
November. 


Compensation practices and local resident perceptions 

We tracked compensation practices back to 2001. The 
overall compensation ratio (51.1%+28.4%) from 2001 to 2014 
fluctuated among schemes, with 72.8% for provincially-funded 


schemes, 85.3% for state-funded schemes, and 93.4% for accounting for annual household income dropped from 14.2% 
insurance-based schemes. The compensation interval was to 4.4% from 2009 to 2014. 
145+33 days during 2005-2013. The proportion of losses 
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Figure 1 Village-scale risk of conflict around Daxueshan Nature Reserve during 2009-2014 
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Figure 2 Grid-scale risk of conflict around Daxueshan Nature Reserve during 2014—2015 
Pie chart shows proportion of losses. The conflicts of Taoshudi (TSD), Baishudi (BSD), Luchang (LC), Hetaoging (HTQ), Kuzhupeng (KZP), Xinzhai (XZ), 
Hongmaoshu (HMS), and Shijiaoyan (SJY) were situated into specific grids. 
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Figure 3 Monthly conflicts of crop depredation caused by 
Asian black bear, rhesus macaque, and Southeast Asian 
sambar around Daxueshan Nature Reserve 
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Figure 4 Monthly conflicts involving goats killed and beehives 
damaged by the Asian black bear around Daxueshan Nature 
Reserve 


We collected 53 questionnaires that comprised 35.1% of 
households involved in conflicts from 2014 to 2015. Results 
showed that although herds were guarded by herders and 
kept in enclosures at night, unexpected Asiatic black bear 
attacks occurred occasionally, with the carcasses of goats 
subsequently having no market value for the farmers. In 
total, 90.6% of respondents were generally satisfied with the 
insurance-based scheme, with 32.1% slightly unsatisfied with 


the delay or unfair assessment and low compensation ratio. 


Furthermore, 60.37% of respondents were willing to pay extra 
insurance premiums for greater compensation. Despite the 
losses, no respondents had a negative attitude regarding 
conservation around DNR. 


DISCUSSION 


Implication to prevent human-wildlife conflict 
Resources to reduce human-wildlife are limited, so it might not 
be necessary to allocate significant effort to less conflict-prone 
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animals or areas. The DNR is inhabited by an assemblage of 
large- and medium-sized mammals, only several of which are 
involved in human-wildlife conflicts. Among those animals, the 
Asiatic black bear was identified as the most damaging and 
tends to be problematic across much of Asia (Escobar et al., 
2015; Lewis et al., 2015; Li et al., 2013). Thus, reduction 
in human-bear conflicts should be addressed first to reduce 
livestock depredation. 


Human-bear conflicts are distinctive events and affected 
by several factors around DNR, highlighting how the black 
bear trade-off between resource extraction and mortality risk 
(Basille et al., 2009; Goswami et al., 2014; Munshi-South 
et al., 2008). Distance to DNR and proportion of cropland 
showed consistent impact on the patterns of livestock and crop 
depredation, both of which were much more frequent from July 
to November, suggesting that livestock and crop depredation 
by the Asiatic black bear were synchronously related. During 
this period, anthropogenic foods, especially maize, are much 
more available than natural foods, thus impacting Asiatic black 
bear foraging (Lewis et al., 2015). The negative correlation of 
distance to DNR revealed that conflicts were partially caused by 
black bears from within the reserve (Karanth et al., 2013). This 
is because suitable habitat for large-sized mammals outside 
the reserve has been replaced by agriculture and infrastructure, 
with the current human population around DNR being 20 times 
higher than the population in the 1950s. The proportion 
of cropland was a major factor reducing the probability of 
the community being affected by livestock depredation. This 
may be because the Asiatic black bear potentially avoided 
predating goats in areas with higher human disturbance, but 
if anthropogenic crops were available, they would risk foraging 
on crops and prey on goats at the same time. Thus, we 
recommend eliminating food crops near the DNR boundary 
as well as livestock grazing at conflict hotspots and offering 
intense guarding to reduce human-bear conflicts. 


Improving compensation practices 


Although monetary compensation receives mixed outcomes 
(Pettigrew et al., 2012), it is still a widespread and 
straightforward approach to satisfy victims of conflict (Dickman 
et al., 2011). Pettigrew et al. (2012) stated that speed, 
transparency, funds, separate responsibilities, involvement of 
experts or trained locals, and clear guidelines are the key 
components of successful schemes. For the insurance-based 
schemes, the speed of compensation and transparency were 
not complained by the respondents around DNR. In addition, a 
third-party insurance company was in charge of assessing and 
compensating losses and the market value of the losses was 
largely compensated. Consequently, the overall attitude toward 
the scheme was positive, rather than the 20% of market value 
compensation for rubber losses in human-elephant conflicts 
and the remoteness and inaccessibility of the Qinghai- Tibet 
Plateau (Chen et al., 2013; Chen et al., 2016b). Moreover, the 
proportion of losses accounting for annual household income 
decreased to one third from 2009 due to the overall increase in 
income. Our results suggest that the insurance-based scheme 


could be replicated at a broader scale with improvement in 
unfair, untimely, and improper assessment. Nevertheless, the 
insurance-based scheme excluded many animals involved in 
human-wildlife conflicts before 2015, such as the wild boar 
and Indian muntjac, despite the wild boar being identified as 
the most conflict-prone animal following the Asian elephant 
(Elephas maximus), accounting for 38.4% of losses in Yunnan 
(Pettigrew et al., 2012). Therefore, it is necessary to enlarge the 
list of animals for which compensation is offered to the victims 
of conflict. 

Currently, the insurance-based scheme is not risk-based 
and does not include resident participation, which has been 
criticized (Chen et al., 2013). It has been proposed that a 
scheme with local participation would ensure more effective 
management of human-wildlife conflict as it would allow 


residents to take greater responsibilities (Chen et al., 2013). 


Our study showed that a high proportion of respondents 
(60.4%) were willing to pay extra insurance premiums to 
receive greater compensation; therefore, schemes with local 
resident participation would improve future compensation. 


CONCLUSIONS 


Understanding the spatiotemporal patterns of human-wildlife 
conflict is imperative for the cost-efficient implementation of 
avoidance and mitigation measures. We presented the 
spatial heterogeneity of conflict at an operational scale with 


spatially-implicit data from conflict compensation schemes. 


Records from mitigation schemes after conflict could facilitate 
the prevention of conflict in the first place. However, studying 
spatial patterns at a finer scale with spatially-implicit data is still 
challenging; thus, we suggest that spatial information should 
be recorded by noting the specific site coordinates and dates 
of conflicts. 
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ABSTRACT 


Lizards are key amniote models for studying 
organ regeneration. During tail regeneration in 
lizards, blastemas contain sparse granulocytes, 
macrophages, and lymphocytes among the prevalent 
mesenchymal cells. Using transmission electron 
microscopy to examine scarring blastemas after third 
and fourth sequential tail amputations, the number 
of granulocytes, macrophages, and lymphocytes 
increased at 3-4 weeks in comparison to the 
first regeneration. An increase in granulocytes and 
agranulocytes also occurred within a week after 


blastema cauterization during the process of scarring. 


Blood at the third and fourth regeneration also showed 
a significant increase in white blood cells compared 
with that under normal conditions and at the first 
regeneration. The extracellular matrix of the scarring 
blastema, especially after cauterization, was denser 
than that in the normal blastema and numerous 
white blood cells and fibroblasts were surrounded 
by electron-pale, fine fibrinoid material mixed with 
variable collagen fibrils. In addition to previous studies, 
the present observations support the hypothesis that 
an increase in inflammation and immune reactions 
determine scarring rather than regeneration. These 
new findings verify that an immune reaction against 
mesenchymal and epidermal cells of the regenerative 
blastema is one of the main causes for the failure of 
organ regeneration in amniotes. 


Keywords: Lizard; Tail; Repetitive regeneration; Blood; 
Immune cells; Ultrastructure 
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INTRODUCTION 


The regenerating tail of lizards is an outstanding and unique 
case of massive organ regeneration in amniotes, i.e., fully 
terrestrial vertebrates (Alibardi, 2010a; Bellairs € Bryant, 1985; 
Fisher et al., 2012; Gilbert et al., 2013). The process of tail 
regeneration in lizards is a uncommon event given the lack of 
any organ regeneration in other amniotes. During regeneration, 
numerous tissues are reformed and combined into a functional 
tail, albeit simplified in comparison to the original. Therefore, 
lizards are critical vertebrates for analyzing successful versus 
unsuccessful tissue and organ regeneration in amniotes, with 
potential consequences for mammalian regeneration (Alibardi, 
2010a, 2014; Gilbert et al., 2013; Lozito & Tuan, 2017). 
Rare cases of anomalous regeneration have been recorded 
in relation to parasitic infection resulting in impaired health 
conditions (Oppliger & Clobert, 1997) or following intense 
damage accompanied by severe inflammation (Baffoni, 1950). 
Tail regeneration can be altered or inhibited under experimental 
manipulations such as cauterization, wounding the regenerating 
blastema, or removing the apical epidermal peg, thus resulting 
in scarring (Alibardi, 2010a, 2013, 2014). Therefore, normal 
regeneration in lizards depends on processes that control 
inflammation but can be interrupted by the above treatments. 
Inflammation is a typical cell and humoral response elicited 
by traumatic mechanical (cuts, amputation), chemical (irritants, 
toxins), or physical (heat, cold, radiation) stimuli to tissues. 
Inflammation limits tissue and organ plasticity in different 
models of vertebrates (Leavitt et al., 2016; Mescher et al., 
2013; Montali, 1988). Invasion of the regenerating tail 
blastema from immune cells, granulocytes, macrophages, 
and lymphocytes is poorly known in lizards, and based 
only on qualitative observations (Alibardi, 1994, 2010b, 2013, 
2016; Alibardi & Sala, 1988). These previous ultrastructural 
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and immunohistochemical studies have indicated that the 
scarring limb undergoes strong invasion from hematogenous 
cells, in particular mononucleate cells (macrophages and 
lymphocytes). 

The importance of clarifying the role of the immune system 
in regenerative plasticity in amniotes was confirmed recently 
based on a transcriptome study comparing the regenerating 
tail blastema to that present in scarring limbs of the same 
lizards (Vitulo et al., 2017). The study indicated that 
inflammation-eliciting genes were unchanged but immune-genes 
were strongly down-regulated in the tail blastema, whereas 
inflammatory genes were up-regulated but immune-genes 
remained unchanged in the limb of the same animals. 

The regenerating blastema therefore appears to function as 
an immune-privileged/suppressed organ, like the developing 
tail or limb, and the presence of immune cells under normal 
tail regeneration is likely limited and does not influence the 
process. Previous electrophoretic and hematological studies 
have indicated that the number of lymphocytes/macrophages 
and the gamma-globulin fraction tend to increase when tail 
regeneration is reduced or inhibited (Alibardi, 2014). However, 
the specific production of anti-regenerative antibodies or 
lymphocytes has not been demonstrated, and therefore the 
increase in the number of immune cells and gamma globulins 
might simply reflect a general state of inflammation. To date, 
no studies have measured the number or type of white blood 
cells, especially macrophages and lymphocytes, invading the 
blastema and circulating in the blood during manipulations 
that stimulate scarring (Alibardi, 2013, 2014). The present 
ultrastructural study aimed to investigate the hypothesis that 
the significant increase in white blood cells, in particular 
macrophages and lymphocytes, likely exerts an immunological 
reaction against blastema cells induced to form scarring 
fibroblasts instead of a mesenchymal mass of cells. 


MATERIALS AND METHODS 


Repetitive tail amputations 

Five adult wall lizards (Podarcis muralis) of both sexes were 
utilized in the present study; they all survived and were 
eventually released in the wild. All experiments followed Italian 


regulations on animal care and handling (art. 5, DL 116/92). 


The lizards were left for 20 min at 4-6 °C to induce hypothermic 
numbness, and the tail was amputated with a razor blade at 
about the second to third proximal. The surface of the stump 
was kept wet and warm with the application of a small drop 
of Ringer fluid, and the extravasating blood from the stump 
was collected with a Pasteur's pipette (30-50 uL) into a small 


0.5-mL Eppendorf vial until no further bleeding was observed. 


A fixative consisting of 0.2 mL of 4% paraformaldehyde solution 
in 0.1 mol/L phosphate buffer at pH 7.4 was rapidly added to the 
blood samples for 5 min. 

The Eppendorf vials were centrifuged for 5 min at 10 000 
r/min to form a pellet, and fixation lasted about 8 h. From 
the same lizards, a 2-mm portion of the normal stump of the 
detached tail was fixed in the same fixative to represent normal 
controls. When the blastema of the first regeneration was 2-3 
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mm in length and coniform in shape, it was cut at the base and 
the blastema with extravasating blood was collected as above 
from the five animals, representing blastema and blood samples 
of the first regeneration. The operation was repeated in the 
following weeks to collect blastema and blood samples of the 
second, third, and fourth regeneration. As controls, three lizards 
were left to regenerate their tails normally. All outgrowths and 
collected blood samples were fixed as above, dehydrated, and 
embedded in Lowcryl resin under UV exposure for 2 d at 4 °C. 


Tail cauterization 

Tail cauterization was performed on three lizards with 
regenerating blastemas 2-3 mm in length. Using a hot plate, 
a laminar spatula was heated for 10 s, then applied three times 
for 1-2 s to the tip of the blastema, scolding it without apparent 
external damage to the skin (i.e., no bleeding or cutting). The 
animals were left for 5-6 d in a cage, with the blastema then 
collected and fixed after growing ceased. The cauterized 
blastemas were cut in half, fixed in 4% paraformaldehyde as 
above and dehydrated in ethanol; one half was then embedded 
in wax (for other studies) and the other half was embedded in 
plastic and transparent resin Lowcryl. 


Preparation for histology and electron microscopy 


Lowcryl-embedded blastemas were sectioned along the 
longitudinal plane using an ultramicrotome and were collected 
on glass slides, stained with 1% toluidine blue, and observed 
under a light microscope for histological analysis. The 
embedded blood samples (pellets) were also sectioned with 
an ultramicrotome. From regions of interest, thin 50—90-nm 
sections were collected on copper grids for transmission 
electron microscopy (TEM) study. The sections were stained 
for 45 min in 1% uranyl acetated water solution and for 5 min 
in lead citrate according to routine methods, then observed 
under a Zeiss C10 electron microscope operating at 60 kv. 
Images were acquired using a digital camera incorporated in 
the electron microscope. 

Two regenerative blastemas (3 mm in length) were prepared 
for the Gomori reaction to reveal acid phosphatases, typical 
markers used to identify lysosome granules in white blood cells 
(Troyer, 1980), under electron microscopy. Tissues were fixed 
at 0—4 °C in 0.25% glutaraldehyde and 4% paraformaldehyde 
in 0.12 mol/L phosphate buffer at pH 7.4. The tissues were 
then cut with a vibratome to collect 50—70-um sections, which 
were incubated for 45 min in 0.5 mol/L acetate buffer at pH 4.8 
and then in Gomori-lead medium for 90 min at 27-30 °C for the 
detection of acid phosphatases (Troyer, 1980). After washing in 
acetate buffer for 30 min, the sections were rinsed in distilled 
water for 5 min, post-fixed in 1% osmium tetroxide for 1 h, and 
then dehydrated and embedded in EPON resin. Thin sections 
obtained using an ultramicrotome were collected on copper grids 
and observed under a Zeiss C10 electron microscope as above. 


Cell counts during successive amputations and after 
cauterization 

Based on the relative percentages of cells in different stages of 
regeneration, two approaches were used to compare the relative 


variations in white blood cells present in blastemas or in blood 
following repetitive amputations of the tail or its cauterization. 
Firstly, the numbers of white blood cells (granulocytes and 
agranulocytes) were counted from whole TEM thin sections 
of regenerative blastemas at the first and fourth regeneration, 
and in cauterized blastemas. The number of detected white 
blood cells in each case was indicated as a percentage out of 
400-500 blastema (mesenchymal) cells counted in blastemas at 
the first and fourth regeneration and one week after cauterization. 
Secondly, the blood samples were examined and the number of 
white blood cells, indicated as a percentage out of 700-900 red 
blood cells, was counted in the blood pellet sections collected 
from lizards in normal, first, third, and fourth regenerations. The 
number of TEM sections (n) for each individual, as well as the 
mean, standard deviation (mean+SD), and statistical significance 
using one-way ANOVA comparing variations in different stages, 
are reported in the text. Statistical significance was set at *: 
P<0.05, **: P<0.01, or ***: P<0.001. 


RESULTS 


Gross and general microscopy observations 

The three lizards allowed to regenerate their tails normally 
produced long and scaled regenerated tails in about 30 d 
following amputation (Figure 1A). Conversely, compared with 
first and second regeneration blastemas, some third and most 
fourth regeneration blastemas appeared as short and hard 
scarring outgrowths. These outgrowths showed irregular scalation 
patterns and were 2-3 mm in length (inset in Figure 1A). 

The normal regenerating tail blastema contained a 
homogenous mass of mesenchymal cells and fibroblasts 
localized underneath a multilayered wound (regenerating) 
epidermis (Figure 1B). Electron microscopy showed that these 
fioroblast-mesenchymal cells possessed a variable number of 
narrow endoplasmic cisternae and euchromatic nuclei with large 
nucleoli. The cells were surrounded by a loose extracellular 
space containing sparse amorphous material (hyaluronic acid) 
and scarce unbanded collagen fibrils (Figure 1C). 

The inner tissues present in the regenerating outgrowths 
at the fourth regeneration displayed less regular connective 
tissue and irregular bundles of extracellular (collagenous) 
fibrils located among fusiform fibroblasts, especially in the 
dermis (Figure 1D). In cauterized blastemas, extracellular fibrils 
or masses of fibrin material were present among fusiform 
fibroblasts and blood vessels (Figure 1E). 


Electron microscopy and count variations at first and 
fourth regenerations 

The normal blastema (first regeneration) contained loose 
connective tissue with many mesenchymal cells and sparse 
macrophages (Figure 1C, F). The Gomori reaction showed 
that macrophages in apical or more proximal regions of 
the blastema contained acid phosphatase positive lysosomes 
(Figures 1G, 2A). In the regenerating outgrowths following 
the fourth regeneration, electron microscopy revealed frequent 
macrophages and granulocytes containing lysosomes in 
various stages of digestion (Figure 1F). In the macrophages, 


lysosomes were acid phosphatase-positive (Figures 1G, 2A), 
and in the granulocytes some granules appeared positive for 
the enzyme as well (Figure 2B). Among the fibroblasts present 
in the outgrowth, numerous collagen fibrils were seen, and the 
extracellular soaces were reduced to form dense connective 
tissues with respect to that of the normal blastema (Figure 2C). 

Furthermore, lymphocytes appeared relatively frequently 
among cells of the blastema in outgrowths at the fourth 
regeneration (Figures 3, 4). Cells recognized as lymphocytes 
were smaller than blastema cells and showed a more 
electron-dense cytoplasm compared with the fibroblasts, with 
variable perinuclear chromatin. Large cytoplasmic regions of 
the lymphocytes contacted the membrane of the blastema 
cells (Figure 3A—C) or occurred through their thin elongation 
(Figure 3D, E). In some cases, the activated lymphocytes 
presented a very irregular surface with typical filopodia, as 
in the macrophages (Figure 4A). In these cases, their exact 
identification was doubtful, and therefore in our quantitative 
counts we considered these cells collectively as agranulocytes 
(see Figure 5). We considered true macrophages only for larger 
cells exceeding 10 um. Other cells with a paler and vesicular 
cytoplasm, but containing large phagosomes or cell remnants, 
were considered as dying phagocytes (Figure 4B). The more 
frequent granulocytes encountered among the connective 
tissues of the outgrowths were represented by basophilic cells. 
These cells were recognized by their numerous and large 
globular granules and by their roundish, non-polymorphic nuclei, 
as observed for neutrophilic granulocytes (Figure 4C). These 
cells were in close contact with the blastema cells, although some 
amorphous material and collagen fibrils were often seen between 
the two cell types (Figure 4D). Small vesicles (0.1—0.2 um in size) 
were also seen in these regions, though their origins remained 
unclear from our static images. 

The number of white blood cells present in first regeneration 
blastemas, out of 400-500 blastema cells, compared to the 
number of white blood cells in fourth regeneration blastemas, 
is presented in Figure 5. Results showed an almost three-fold 
increase in white blood cells at the fourth regeneration (mean 
15.9), especially agranulocytes, compared to that at the first 
regeneration (mean 5.0). This value was significant according 
to the ANOVA test (Table 1, Figure 5). 


Electron microscopy and count variations in cauterized 
outgrowths 

After cauterization, the blastema stopped its growth. Among 
blastema cells (fibroblasts), a reduced extracellular matrix was 
present in comparison to the normal blastema (Figure 6A, B). 
White blood cells were more frequently encountered than in 
the normal blastema at the first regeneration and appeared as 
agranulocytes (macrophages-lymphocytes) and granulocytes 
with large granules, likely basophils. Quantification indicated a 
three-fold increase (mean 15.1, see Figure 5) in white hematic 
cells in the blastema at 6 d post-cauterization compared with 
the number of white cells present in the blastemas in the 
first regeneration (mean 5.0, see Figure 5). This value was 
significant according to the ANOVA test (Table 1, Figure 5). 
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Figure 1 Gross aspects of regenerated tail versus scarring tail outgrowth (A) and microscopy images of areas of the regenerating 
tail (B—G) 

A: Largely scaled, regenerated tail at 30 d post-amputation. Arrowhead indicates transition point with stump of original tail. Scale bar: 1 mm. Inset shows a 
tail outgrowth at fourth regeneration (arrowhead on amputation region). Scale bar: 0.5 mm. B: Histology of a first regeneration blastema. Scale bar: 10 um. 
C: Electron microscopy view of first regeneration blastema cells. Arrows indicate collagenous fibrils with amorphous matrix. Scale bar: 0.5 um. D: Histology 
of a fourth regeneration blastema. Scale bar: 10 um. E: Histology of a blastema at 6 d post-cauterization, showing large masses of pale extracellular material 
(arrows) underneath wound epidermis. Scale bar: 10 um. F: Macrophage with numerous secondary lysosomes (arrowheads) and filopodia (arrows) within a 
first regeneration blastema. Scale bar: 1 um. G: Macrophage within the apical region of a first regeneration blastema with acid phosphatase positive lysosomes 


(arrows). Scale bar: 1 um. de: dermis; ex: extracellular matrix; n: nucleus; w: wound epidermis. 
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Figure 2 Electron microscopy view of white blood cells in a first regeneration blastema (A, B) and a fourth regeneration 


fibroblast (C) 
A: Macrophages with acid phosphatase positive granules (arrows) in more proximal area of the blastema. Scale bar: 0.5 um. B: Granulocyte with some acid 


phosphatase positive granules (arrows). Scale bar: 0.25 um. C: Fibroblast in contact with large bundle of collagen fibrils. Scale bar: 1 um. col: collagen fibrils; 


ex: extracellular matrix; n: nucleus; ns: non-specific granules. 
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Figure 3 White blood cells present in a fourth regeneration blastema 

A: Mesenchymal fibroblast close to darker agranulocyte (arrow). Scale bar: 1 um. B: Large appositional area (arrowheads) between lymphocyte (darker) and 
blastema cells (paler fibroblast). Scale bar: 0.25 um. C: Contact area with filopodia (arrows) localized between a lymphocyte and fibroblast. Scale bar: 0.25 um. 
D: Macrophage filopod contacting (arrow) mesenchymal cell. Scale bar: 0.25 um. E: Filopod from a macrophage contacting or invading (arrows) a blastema 


cell. Arrowheads point to seemingly extracellular vesicles. Scale bar: 0.25 um. ex: extracellular matrix; ly: lymphocyte; ma: macrophage; n: nucleus. 


A pale and finely granular material, resembling a filaments (Odland & Ross, 1968), were also embedded in the 
plasma-fibrinous exudate, was present in the extracellular exudate. These fibrils were irregularly distributed within the 
matrix (Figure 6C—E). Denser coarse filaments or fibrils of extracellular matrix. 


30-40 nm and of unknown nature but resembling fibrin 
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Figure 4 Electron microscopy view of white blood cells within a fourth regeneration blastema 


A: Agranular electron-dense cell with filopodia (arrows) and cytoplasmic dense bodies (arrowheads, likely lysosomes). Scale bar: 1 um. B: Pale cell with two 


large phagosomes (arrowheads) near a cell with filopodia (arrows). Scale bar: 1 um. C: Basophil granulocyte (arrows show large dense granules) detected 


among blastema cells. Scale bar: 1 um. D: Narrow extracellular soace present between a basophil and blastema cell, where collagen fibrils and amorphous 


material are present (arrows). Arrowheads point to apparent extracellular vesicles. Scale bar: 0.25 um. bl: blastema (mesenchymal) cell; col: collagen fibrils; 


ex: extracellular space; n: nucleus. 


Electron microscopy of blood cells in normal, first, third, 
and fourth regeneration 

The blood cells derived from pellets of normal blood samples 
and after the first, third, and fourth regeneration demonstrated 
a significant 3—4-fold increase in white blood cells in general 
(1.2% in normal blood versus 5.1% in blood at the fourth 


regeneration) (Table 2, Figure 7). The value was not significant 
between normal and first regeneration or normal and second 
regeneration but was significant between first and second 
regeneration versus fourth regeneration and highly significant 
between normal and fourth regeneration (Figure 7). 
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2 - White blood cells in both pellets and blood vessels 
in the blastema were mainly represented by granulocytes 
(especially basophils), though macrophages and lymphocytes 
were also seen (Figure 8A). Numerous cells were likely 
lymphoblasts, characterized by a round shape and pale 
cytoplasm, and containing sparse free ribosomes and a 
nucleus featuring a more diffuse (unpacked) chromatin than 
that of small lymphocytes, macrophages, or granulocytes 
(Figure 8B—E). Completely differentiated plasma-cells were 
only found occasionally in blood pellets sections and within 
vessels in the blastema or in scarring outgrowths (Figure 8F). 
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white blood cells/blastema cells in first regeneration (1stR) 
compared to fourth regeneration (4thR) or cauterized (Caut) 
**: P<0.01; ***: P<0.001. 
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Figure 6 Electron microscopy view of blastema 6 d after cauterization 
A: Basophil granulocytes (arrow) and macrophages (arrowhead) present among blastema cells. Scale bar: 2 um. B: Electron-pale extracellular matrix present 


among blastema cells. Scale bar: 1 um. C: Dark lymphocyte within the electron-pale extracellular exudate. Scale bar: 1 um. D: Dense amorphous material 
present within the extracellular pale exudate. Scale bar: 0.25 um. E: Denser fibrillar material (arrows) embedded in paler exudate between two dark cells 


(arrowheads). Scale bar: 0.2 um. bl: blastema cells; dsm: dense extracellular and largely amorphous material; ex: extracellular space; n: nucleus. 
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DISCUSSION 


This ultrastructural study demonstrated a significant increase 
in white blood cells, about three times higher than that 
observed during the first regeneration, in the cauterized and 
fourth repetitive regeneration blastemas, confirming previous 
studies (Alibardi, 2010a, 2010b, 2013, 2014). In particular, 
mononucleate macrophages or activated lymphocytes were 
predominant over other blood cell types (Figure 5), whereas 
granulocytes were more abundant than agranulated cells in 
blood (Figure 7). Wound exudate, which contained fibrin, was 
abundant in the cauterized blastemas, likely due to plasma 
extravasation following cauterization. The removal of the 
exudate likely generated intense inflammation in the blastema, 
leading to scarring (Alibardi, 2013). 


Table 1 Percentages of white blood cells in first regenerative 
blastema compared with those present in fourth regenerative 
blastema and 6 d after blastema cauterization 


1streg 4th reg Cauterized 
n 5 4 3 
G (%) 45.2 28.8 40.2 
A (%) 54.8 71.2 59.8 


WBC (%) 5.0+1.5 15.941.6 15.1+1.8 


n: Number of samples. G: Granulocytes (heterophils-neutrophils, 
basophils, and eosinophils). A: Agranulocytes (monocytes and 
lymphocytes). WBC: White blood cells. 


Table 2 Percentages of white blood cells in the blood of normal 
lizards and in first, third, and fourth regeneration 


Normal 1st reg 3th reg 4th reg 





n 3 4 3 3 
G (%) 56.7 56.2 58.0 54.9 
A (%) 43.3 43.8 42.0 45.1 


WBC (%) 1.2+0.45 2.140.41 2.540.387 5.1+0.53 


n: Number of samples. G: Granulocytes (heterophils-neutrophils, 
basophils, and eosinophils). A: Agranulocytes (monocytes and 
lymphocytes). WBC: White blood cells. 


Most lymphocytes contained sparse ribosomes and lacked 
endoplasmic reticulum vesicles, suggesting they were T-cells 
(Janossy et al., 1972). However, due to some uncertainty in the 
discrimination between macrophages and activated lymphocytes 
or between activated basophil and eosinophil granulocytes, the 
present quantitative data on white cell tyoes should be considered 
with some caution. In the average sections of the blastema 
or blood, cells appeared in their active, often motile condition, 
and their shape and surface were distorted and not optimal for 
unbiased morphological identification. Furthermore, the plane 
of the randomly cut sections may also have intercepted cells 
in a cytoplasmic or nuclear region that were tangential to the 
main cellular area, making identification more problematic. Some 
studies have also indicated that lymphocytes in fish, amphibians, 
and some reptiles possess active phagocytosis, making the 
distinction with macrophages even more difficult (Li et al., 2006; 


Zimmerman et al., 2010). Despite these shortcomings, the present 
study clearly showed that significant increases in the number of 
granulocytes, lymphocytes, and particularly macrophages in the 
blastema and blood were associated with scarring. 
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Figure 7 Statistical values (ANOVA test) and significance of 


variations in white blood cells/erythrocytes in normal blood 
(NB), first regeneration (1stR), third regeneration (3rdR), and 
fourth regeneration (4thR) 

**- P<0.01; ***: P<0.001. 


The present study indicated that an inflammatory reaction 
involving an average of 15% of white blood cells invading the 
blastema was strong enough to limit or impair tail regeneration 
in lizards. Therefore, it is not surprising that the healing 
limb of Podarcis muralis at 11-12 d of regeneration after 
amputation, where inflammation includes over 50% of white 
blood cells in comparison to the total cell types present 
in a forming blastema (Alibardi, 2010b, 2016; 50.5%+2.4% 
in three quantified cases), gives rise to a scar. The 
main proportion of immune cells, about 70% of the total, 
consisted of macrophages/lymphocytes of difficult specific 
identification. The increase in basophils/eosinophils and 
macrophages/lymphocytes detected in the present study in 
scarring tails indicated that these cells were likely associated 
with the formation of scarring connective tissue. However, 
direct stimulation of myofibroblasts and fibrocytes in lizards 
from these cells was not shown. The differentiation of the 
latter cells can give rise to a scarring blastema during the fourth 
regeneration in a cauterized tail stump (Alibardi, 2010a, 2013) 
or in the cauterized blastema (present study). 

Numerous studies on mammalian and non-mammalian 
vertebrates have shown that inflammation, especially chronic 
inflammation, is detrimental to regeneration (Ferguson & 
O’Kane, 2004; Leavitt et al., 2016; Mescher et al., 2013). 
However, the type of injury can evocate different types of 
inflammatory cells, in particular macrophages, some of which 
(healing or M2) actually favor regeneration in amphibians 
(Godwin & Rosenthal, 2014), fish (Petrie et al., 2014), and 
mammals (Hesketh et al., 2017; Simkin et al., 2017). The best 
regenerators among vertebrates, urodele amphibians, and also 
many fish, possess low inflammation and immune responses 
after injury, Supporting the hypothesis that immune efficiency is 
the main obstacle to organ regeneration in amniotes. Another 
clear indication of the negative effect of the immune system 
on organ regeneration is present in anuran amphibians, which 
are more terrestrially-adapted than many urodeles (Mescher et 
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al., 2013, 2017). After metamorphosis, the immune response 
is enhanced, and the regeneration capability is lost in frogs, 
toads, and some terrestrial salamanders (Mescher et al., 2013, 
2017; Scadding, 1977, 1981). Evolution and potentiation of 
the adaptive immune system are correlated with the almost 


complete loss of organ regeneration in endothermic amniotes. 


The strong inflammation elicited in both avian and mammalian 
organs after injury leads to rapid resolution and variable 
degree of scarring (Ferguson & O'Kane, 2004; Hesketh et al., 
2017). It is likely that after limb amputation or cauterization in 
lizards, traumatic events that induce considerable inflammation, 


the type of macrophages and lymphocytes activated are 
pro-inflammatory (M1), leading to scar formation. 

The present quantitative ultrastructural study supports the 
hypothesis that the immune system is likely the main cause 
of failure of tissue and organ regeneration in lizards, and in 
amniotes in general (Alibardi, 2014). Future follow-up studies 
on this topic are required to demonstrate whether inflammatory 
cells respond to general stimulation after intense and persistent 
tissue damage or if macrophages and lymphocytes can 
specifically attach to mesenchymal, ependymal, and epidermal 
cells of the regenerating blastema. 





Figure 8 Electron microscopy view of blood cells present in blood vessels within a blastema (A, D—F) or isolated in blood pellets 


(B, C) at the fourth regeneration 


A: In this large capillary, most lymphocytes (arrowheads), lymphoblasts (arrows), and monocytes (double arrows) are seen. Scale bar: 2 um. B: Pale lymphoblast 


detected among erythrocytes. Scale bar: 1 um. C: Lymphoblast present among erythrocytes. Scale bar: 1 um. D: Polymorphonucleate granulocyte (arrowheads 


indicate three lobes of the nucleus). Scale bar: 1 um. E: Mastocyte (basophil?) with filopodia (arrows). Scale bar: 1 um. F: Plasma cell within blastema. Scale 


bar: 2 um. col: collagen bundle; end: endothelial cell; ex: extracellular matrix; ery: erythrocyte; n: nucleus; ns: non-specific granules; pla: plasma-cell. 
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ABSTRACT 
The Bama Xiang pig (BMX) is a famous early-maturing 


Chinese indigenous breed with a two-end black coat. 


To uncover the genetic basis of the BMX phenotype, 
we conducted comparative genomic analyses between 
BMX and East Asian wild boars and Laiwu pigs, 
respectively. Genes under positive selection were 
enriched in pathways associated with gonadal hormone 
and melanin synthesis, consistent with the phenotypic 
changes observed during development in BMX pigs. We 
also performed differentially expressed gene analysis 
based on RNA-seg data from pituitary tissues of BMX 
and Large White pigs. The CTTNBP2NL, FRS2, 
KANK4, and KATNAL7 genes were under selection 
and exhibited expressional changes in the pituitary 
tissue, which may affect BMX pig puberty. Our study 
demonstrated the positive selection of early maturity 
in the development of BMX pigs and advances our 
knowledge on the role of regulatory elements in puberty 
evolution in pigs. 


Keywords: Puberty; Bama Xiang pig; Pituitary; 
Differentially expressed genes 


INTRODUCTION 


Pigs are important protein food source and globally distributed 
animal domesticated from Eurasian wild boars 9 000 years ago 
(Larson et al., 2005). During their long domestication history, 
pigs have exhibited abundant phenotypic variation in coat color, 
reproduction, and growth. This phenotypic diversity provides 
valuable animal models for genetic studies (Ciobanu et al., 
2001; Rubin et al., 2012). 

The Bama Xiang pig (BMX) from Guangxi Zhuang 
Autonomous Region is a well-known variety noted for its early 
maturity and two-end black pigmentation. The female BMX 
pig matures at about three months, in sharp contrast to the 
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five months observed in most other East Asian domestic pigs 
(China National Commission of Animal Genetic Resources, 
2011). Thus, the BMX pig provides an important model to 
study early maturity and pigmentation evolution during pig 
domestication. 

In this study, we performed genomic screening of signatures 
of positive selection in the BMX pig in comparison with those 
of the East Asian wild boar (EAW) and Laiwu (LWU) pigs. The 
LWU pig has a normal maturity process (at six months) and 
black coat color (China National Commission of Animal Genetic 
Resources, 2011), but shares the same domestication origin as 
that of the BMX pig from East Asian wild boars (Wu et al., 2007). 
In combination with transcriptomic analysis, we identified genetic 
variations that were putatively associated with the maturity and 
pigment phenotype in the BMX pig. The primary goal of this 
study was to explore the (1) genomic evolutionary pattern shaped 
by development of the BMX pig, and (2) genetic variations 
associated with selection for early maturity and coat color in the 
BMX pig. 


MATERIALS AND METHODS 


Sample collection and resequencing of genomic DNA 

We downloaded the genome resequencing data of six BMX, 
six LWU, and 13 EAW pigs from the NCBI database 
(Supplementary Table S1), as reported from earlier studies 
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(Ai et al., 2015; Groenenet al., 2012; Li et al., 2013; Rubin 
et al., 2012). Sequence reads were filtered by removing 
adaptors and low-quality bases using Trimmomatic (Bolger et 
al., 2014). Genomic reads were aligned to the Sus scrofa 10.2 
genome (Groenenet al., 2012) with the BWA program (Li & 
Durbin, 2010). Alignment was performed with a maximum of 
six mismatches allowed for each 100-bp fragment. 


Single nucleotide polymorphism (SNP) calling and annotation 
SNP calling was performed with the DNA sequencing data 
of the 25 pig genomes using HaplotypeCaller in the GATK 
toolkit (McKenna et al., 2010). Only polymorphic SNPs were 
included in this analysis. VariantFiltration and SelectVariants 
were used to filter SNPs by quality score (QUAL <20) and depth 
of coverage (DP <10). SNPs were annotated with the GTF file 


downloaded from the NCBI website (https://www.ncbi.nlm.nih. 


gov/). For annotation, SNPs were classified by their locations 
in protein-coding (Synonymous or non-synonymous), intronic, 
intergenic, and untranslated region (UTR) sequences. 


Selective sweep detection 

Population fixation index (Fst) (Weir 8 Cockerham, 1984), cross 
population extended haplotype homozygosity (XPEHH) (Sabeti 
et al., 2007), and nucleotide diversity (x) (Tajima, 1983) were 


used to search for signals of selection in the BMX pig genome. 


We analyzed selective sweeps in 10-kb non-overlapping sliding 
windows. Windows containing <5 SNPs were removed from the 
analysis. 

The Fs; values between BMX and other pigs (LWU and 
EAW) were calculated using VCFtools (Danecek et al., 2011). 

Haplotypes were inferred using Beagle software (Browning 
& Browning, 2007). Comparing the BMX pigs with LWU and 
EAW pigs, XPEHH was calculated to examine haplotypes that 
showed low levels of linkage decay in the BMX populations 
(Coop et al., 2009; Pickrell et al., 2009; Pritchard et al., 2010). 

The nucleotide diversities of BMX (gmx), EAW (meaw), 
and LWU (mıwu) pigs were then calculated separately using 
VCFtools (Danecek et al., 2011). To identify regions that 
showed low nucleotide diversity in the BMX pigs compared to 
the LWU and EAW pigs, we calculated meaw divided by tpux 
and TLWU divided by TBMX> respectively. 

Genomic windows showing high Fsr, low nucleotide diversity, 
and high XPEHH values were identified as having signatures 
of selection in BMX pigs. Candidate selective sweeps were 
chosen from genomic windows identified simultaneously by all 
three statistics as having values above the top 10% threshold in 
their empirical distributions. 


RNA sample collection, sequencing, and identification of 
differentially expressed genes 

Pituitary tissues were collected from three BMX pigs and three 
Large White pigs for RNA sequencing (RNA-seq). These pigs 
were females at 85 days after birth. The RNA-seq libraries 
with an inserted size of 250 bp were prepared using the 
Illumina standard RNA-seq library preparation pipeline and 
sequenced on the Illumina Hiseg 2000 platform, with 150-bp 
paired-end reads generated. In total, the transcriptomes of 


six pituitaries (GSA number: CRA000876) were used for 
comparative analysis. The RNA-seg clean reads were mapped 
onto the Sus scrofa 10.2 genome using Tophat2 (Kim et al., 
2013). Anew merged GTF annotation file was generated using 
Cufflinks and Cuffmerge (Trapnell et al., 2010). 

Differentially expressed genes (DEGs) were identified by 
DESeq2 (Love et al., 2014). As DESeq2 required a raw read 
count table as input, we obtained the raw read count using 
HTseq (Anders et al., 2014). An FDR of <0.05 was used as the 
cutoff for DEGs. The ClueGO plugin in Cytoscape (Mlecnik et 
al., 2018) was used for gene ontology (GO) enrichment analysis. 
GO terms with a P-value of <0.05 were defined as enriched. 


RESULTS 


Selective signatures in BMX pig genomes 

The whole genomes of six BMX, 13 EAW, and six LWU 
pigs were used to identify genomic variants underlying the 
development of BMX pigs. Genomic reads were quality filtered 
and mapped to the pig reference genome (Sus scrofa 10.2) 
(see Materials and Methods). A total of 43.5 million SNPs 
were identified in the 25 individuals. To identify signatures of 
selection in BMX pig genomes, we analyzed the Fs; (Weir 
& Cockerham, 1984), x (Tajima, 1983), and XPEHH (Sabeti 
et al., 2007) metrics along chromosomes in comparison to 
those of EAW and LWU pigs (Figure 1). Selective signatures 
were screened in 10-kb non-overlapping sliding windows and 
adjacent selective windows were combined. As recent artificial 
selection produces long stretches of haplotypes in a population 
(Amaral et al., 2008), a clustering of selective windows with 
intervals of less than nine windows was chosen as an indicator 
of a putatively selective sweep. Sweep regions greater than 
30 kb in length were included in the following analyses. This 
procedure resulted in the identification of 257 genes within 
339 putative selective sweeps in the BMX pig genomes 
(Supplementary Table S2). 

Functional enrichment analysis of the gene set revealed 
significant overrepresentation of 59 terms (Table 1). 
Interestingly, several biological processes related to sexual 
maturity were involved, including “pituitary gland development” 
(BMP4, FGF10, and TBX79), “endocrine system development” 
(BMP4, CRH, FGF10, NEUROD1, and TBX19), and “negative 
regulation of reproductive process” (BMP4, CSN2, and ZP4). 
The pituitary is the most important endocrine gland involved 
in growth and reproduction regulation. In addition, five genes 
(FGF10, FGF2, MITF, PDGFC, and PDGFD) were involved 
in “Melanoma”, which might be associated with the two-end 
black pigmentation observed in BMX pigs, which differs from 
most Chinese domestic pigs with black coats (China National 
Commission of Animal Genetic Resources, 2011). 

Although many missense mutations are associated with the 
domestication of pigs (Fang et al., 2009; Li et al., 2010; Ren 
et al., 2011; Rubin et al., 2012), whether regulatory variants 
play a role and their effects on the domestication processes 
compared to coding variants remain unknown. To address 
this issue, we compared the evolutionary patterns for SNPs 
identified from the selective regions in the BMX pig genomes. A 
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total of 830 949 SNPs were identified from 339 sweeps. Among 
these, 3 668 SNPs were located in coding sequences and 
827 542 SNPs were located in noncoding sequences. We 
first focused on protein coding region variations located in 
the 339 sweeps. We examined highly differentiated missense 
mutations and analyzed their presence in DNA sequences 
encoding conserved protein sequences. We concentrated on 
missense SNPs with Fs; values up to the chromosomal top 
0.1% threshold and compared BMX pigs to LWU and EAW 
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pigs, respectively. We identified a dataset consisting of 11 
missense SNPs from 10 genes (Supplementary Table S3). 
Alignment of the 10 protein sequences showed that the flanking 
sequences at the residues encoded by the 11 SNPs displayed 
low level cross-species conservation, implying weak functional 
constraints on these protein-coding variants. This indicated 
that regulatory elements rather than missense variants may 
have greater effects on BMX pig development. 















Figure 1 Whole genome selective signatures of BMX pigs 


X-axis shows the positions of the windows along the 1-18 autosome chromosomes. Different chromosomes are displayed in different colors. Y-axis shows 


the Fst (A), XPEHH (B), teawemx (C), and tLwu/pmx (D) values. Fst and XPEHH values were calculated between BMX and (LWU+EAW) pigs. teawemx and 


TLwu/emx Were obtained from meaw and mwy divided by TBmx, respectively. 


Comparative transcriptomic analysis of pituitary tissues 

We conducted transcriptomic profiling to examine gene 
expression changes due to selection in the BMX genomes as 
regulatory variants exhibited a profound role during BMX pig 
development. We chose pituitary tissue as the pituitary gland 
secretes reproduction-related hormones involved in maturity 
regulation. In our experimental design, we collected pituitary 
tissues from three BMX and three Large White pigs at 85 
days after birth. At this stage, BMX pigs are close to maturity, 
whereas Large White pigs do not reach maturity until eight 
months of age (China National Commission of Animal Genetic 
Resources, 2011). The transcriptomes from the pituitary tissues 
were sequenced, from which we obtained 35 million mapped 
reads, on average, for each individual. We compared the 
gene expression levels in the pituitary and found 1 600 DEGs 
between BMX and Large White pigs. Although Large White 
pigs were domesticated from European wild boars (Larson et 
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al., 2005;Rubin et al., 2012), the DEGs between BMX and 
Large White pigs provided valuable information on expressional 
changes in BMX pigs if they also exhibited selective signatures 
in the BMX genome. The DEGs were enriched in the 
following biological processes: “Estrogen signaling pathway” 
(ADCY6, ADCY8, CREB3L1, CREB5, FOS, GABBR2, GNAS, 
GRM1, KCNJ9, MMP9, PIK3CG, PIK3R5, and PRKACA), 
“Ovarian steroidogenesis” (ADCY6, ADCY8, ALOX5, GNAS, 
LDLR, LHCGR, and PRKACA), “Oocyte meiosis” ( ADCY6, 
ADCY8, CAMK2B, CAMK2D, CCNE1, CDC20, CPEB2, ESPL1, 
MAD2L2, MOS, PRKACA, and YWHAZ), “Developmental 
maturation” ( C12H17orf47, CDC20, EDNRB, FAM101B, LTF, 
PLP1, PRKACA, RBPJ, TAL1, WASH1, and WNT10B), and 
“Progesterone-mediated oocyte maturation” ( ADC Y6, ADC Y8, 
CPEB2, MAD2L2, MOS, PIK3CG, PIK3R5, and PRKACA) 
(Supplementary Table S4). These pathways all participate in 
sexual maturity. 


Table 1 Functional enrichment terms of positively selected genes in BMX pigs. 





GO ID GO term P-value’ 
GO:0004080 Neuroactive ligand-receptor interaction 24.0E-15 
GO:0060449 Bud elongation involved in lung branching 100.0E-6 
GO:0032740 Positive regulation of interleukin-17 production 140.0E-6 
GO:0060393 Regulation of pathway-restricted SMAD protein phosphorylation 230.0E-6 
GO:0060389 Pathway-restricted SMAD protein phosphorylation 250.0E-6 
GO:0032946 Positive regulation of mononuclear cell proliferation 310.0E-6 
GO:0050671 Positive regulation of lymphocyte proliferation 310.0E-6 
GO:0070665 Positive regulation of leukocyte proliferation 370.0E-6 
GO:0048636 Positive regulation of muscle organ development 400.0E-6 
GO:0045844 Positive regulation of striated muscle tissue development 400.0E-6 
GO:1901863 positive regulation of muscle tissue development 450.0E-6 
GO:0005412 Arrhythmogenic right ventricular cardiomyopathy (ARVC) 490.0E-6 
GO:0005218 Melanoma 560.0E-6 
GO:0045740 Positive regulation of DNA replication 600.0E-6 
GO:0009798 Axis specification 600.0E-6 
GO:0090090 Negative regulation of canonical Wnt signaling pathway 600.0E-6 
GO:0032660 Regulation of interleukin-17 production 660.0E-6 
GO:0051155 Positive regulation of striated muscle cell differentiation 660.0E-6 
GO:0032620 Interleukin-17 production 780.0E-6 
GO:0042506 Tyrosine phosphorylation of Stat5 protein 780.0E-6 
GO:0051145 Smooth muscle cell differentiation 800.0E-6 
GO:0061037 Negative regulation of cartilage development 920.0E-6 
GO:0048546 Digestive tract morphogenesis 950.0E-6 
GO:0005410 Hypertrophic cardiomyopathy (HCM) 990.0E-6 
GO:0051149 Positive regulation of muscle cell differentiation 1.2E-3 
GO:0035270 Endocrine system development 1.3E-3 
GO:0055025 Positive regulation of cardiac muscle tissue development 1.3E-3 
GO:0042310 Vasoconstriction 1.5E-3 
GO:0060602 Branch elongation of an epithelium 1.5E-3 
GO:0001 709 Cell fate determination 2.0E-3 
GO:0021536 Diencephalon development 1.8E-3 
GO:0021983 Pituitary gland development 2.2E-3 
GO:0042307 Positive regulation of protein import into nucleus 2.3E-3 
GO:0000561 Glycerolipid metabolism 2.4E-3 
GO:0004730 Long-term depression 2.6E-3 
GO:0003338 Metanephros morphogenesis 2.4E-3 
GO:0030199 Collagen fibril organization 2.4E-3 
GO:0060441 Epithelial tube branching involved in lung morphogenesis 3.0E-3 
GO:0030219 Megakaryocyte differentiation 3.3E-3 
GO:0042383 Sarcolemma 3.3E-3 
GO:0009799 Specification of symmetry 3.3E-3 
GO:0051153 Regulation of striated muscle cell differentiation 3.3E-3 
GO:0010092 Specification of organ identity 3.3E-3 
GO:1900182 Positive regulation of protein localization to nucleus 3.7E-3 
GO:0005033 Nicotine addiction 4.0E-3 
GO:0004140 Regulation of autophagy 5.9E-3 
GO:2000242 Negative regulation of reproductive process 5.9E-3 
GO:0070664 Negative regulation of leukocyte proliferation 5.1E-3 
GO:0003401 Axis elongation 5.1E-3 
GO:0009948 Anterior/posterior axis specification 5.5E-3 
GO:0055024 Regulation of cardiac muscle tissue development 5.1E-3 
GO:0004940 Type | diabetes mellitus 7.4E-3 
GO:0048259 Regulation of receptor-mediated endocytosis 6.9E-3 
GO:0007193 Adenylate cyclase-inhibiting G-protein coupled receptor signaling pathway 11.0E-3 
GO:0009880 Embryonic pattern specification 11.0E-3 
GO:0005160 Transforming growth factor beta receptor binding 9.5E-3 
GO:0031623 Receptor internalization 11.0E-3 
GO:0048286 Lung alveolus development 9.5E-3 
GO:0010862 Positive regulation of pathway-restricted SMAD protein phosphorylation 10.0E-3 


": P-value obtained from ClueGO plugin in Cytoscape using default paraments. Only terms with P<0.05 are listed. 
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As both artificially selected genes and DEGs showed 
functional enrichment in maturity regulation, it is possible 
that differential gene expression was favored by the 
artificial selection of variants from some sweep loci. We 
therefore identified artificially selected genes with altered gene 
expression. In our analysis, a total of 22 artificially selected 
genes exhibited differential gene expression in the pituitary 
(P<0.05; FDR<0.05; Table 2). Interestingly, four genes 
(CTTNBP2NL, FRS2, KANK4, and KATNAL 1) were related to 
reproduction (Figure 2). Compared to that in the Large White 
pigs, the expression of the CTTNBP2NL gene was increased in 
BMX pigs. CT TNBP2NL has been reported to affect conception 
rates in bovines (Sugimoto et al., 2013). FRS2 encodes a 





© 
© 
Wy 
© 
© 
© © 
g T 2 + 
= € 
2 = 
9 S 2 
a o + 
S o [q] 
o Sa D 
= 00 — 
a) O o 
O oO ra) 
A> > m 
E E 
z E o 
o e o Se 
ZE a A 
© 
N 
N 
© 
© © 
+ © 
N 
BMX LargeW BMX LargeW 
CTTNBPINL FRS2 


value, black horizontal line in box is median, and lower whisker is minimum value. 


DISCUSSION 


We conducted a comparative genomic analysis of BMX pigs 
with domestic pigs and wild boars from East Asia, and revealed 
selective signatures potentially associated with selection for 
early maturity and the two-end black phenotype. During 
domestication, selection is based on those functions and traits 
that are favored in pig breeds. As described earlier, selective 
signatures were detected for the specific phenotype of the BMX 
pig. This study sheds light on the evolution of puberty in pigs 
and provides important information and candidate genes for pig 
breeding. 


Puberty is an important trait related to economic output, 
with early maturity resulting in early slow growth. For 
instance, Western commercial breeds are characterized with 
late maturity for enlarged body size (Rubin et al., 2012). The 
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docking protein bound to FGF receptors, with a knockout study 
in mice indicating its involvement in epididymal development 
and sperm maturation (Xu et al., 2013). In the pituitary tissue, 
we found that the expression level of FRS2 in BMX pigs was 
about twice that found in Large White pigs. In addition, KANK4 
was found to be differentially expressed in the pituitary tissues 
of BMX and Large White pigs. Although few studies have 
reported on the function of KANK4, it is located in a significant 
pig quantitative trait locus (QTL) for age at puberty (Bidanel 
et al., 2008). KATNAL1, which encodes Katanin Catalytic 
subunit A1 like 1, has been implicated in both meiosis and 
spermiogenesis (Smith et al., 2012). 
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Figure 2 Normalized read counts of significantly differentially expressed genes CTTNBP2NL (A), FRS2 (B), KANK4 (C), and 
KATNAL 1 (D) in the pituitary 


For each subgraph, red and green boxes represent gene expression in the pituitaries of BMX and Large White pigs, respectively. Upper whisker is maximum 


present study revealed four loci (CT TNBP2NL, FRS2, KANK4, 
and KATNAL1) showing strong signatures of selection and 
functional association with puberty. In addition, the four 
genes exhibited differential expression between early and late 
maturity in pig breeds in the pituitary. This result indicates that 
regulatory elements may have considerable effect on BMX pig 
development. 


Coat color is a trait that is usually under strong selection 
in different breeds. In previous studies, several genes have 
been shown to be associated with coat color variants, including 
MC1R, KIT, and EDNRB (Kijas et al., 1998; Wiener & Wilkinson, 
2011; Wilkinson et al., 2013). However, we found no selective 
signature of MC1R. The reason for this may be that selection of 
MC1R occurred during early domestication, and most Chinese 
indigenous pig breeds share the same haplotype (Li et al., 
2010). Earlier research has indicated that mutations in MITF, 


PAX3, EDNRB, and KIT together may explain a large proportion 


of white spotting phenotypes in horses (Hauswirth et al., 2012). 


In BMX pigs, the two-end black phenotype is a complex trait 
resulting in the diversification of white areas on the body. We 
hypothesize that a correlation exists between coat color variants 
and different combinations of candidate genes. 


Table 2 Overlap between positively selected genes and 
differentially expressed genes from pituitary tissue in 
BMX pigs. 





Gene Normalized count fold change (BMX/LW) P-value’ 
MGAT5B 0.08 7.57E-05 
BHLHE22 2.57 0.000 538 
FAM174B 0.43 0.000 589 
SFT2D2 0.28 0.001 381 

TRHDE 2.13 0.005 172 
COL23A1 0.12 0.006 991 
FRMD4B 2.45 0.007 354 
CDKN2B 0.14 0.007 441 

FRS2 2.09 0.008 626 

ATG4C 2.16 0.010 881 

PELI2 1.90 0.012 313 

FEZ2 2.12 0.018 321 
LGALS12 0.13 0.020 835 

TYRP1 0.13 0.025 025 
KCND3 1.81 0.030 48 
FNDC1 0.22 0.030 808 

ATL3 0.46 0.033 021 

KANK4 0.34 0.036 521 
KATNAL 1 1.99 0.037 825 

CTTNBP2NL 2.27 0.039 048 

ANOS 2.50 0.044 986 

DYNC112 1.67 0.047 03 


": P-values are corrected by Benjamini-Hochberg FDR<0.05. 


One limitation of this study was the identification of the DEG 
direction in RNA-seq analysis. As both BMX and Large White 
pigs are domesticated, it was difficult to identify the direction of 
gene expression without the use of an outgroup. We analyzed 
the overlap between DEGs (between BMX and Large White 
pigs) and the artificially selected genes in BMX pigs (compared 
to LWU and EAW pigs). The shared genes were more likely 
to be changed in BMX pigs. We did not focus on DEGs that 
did not exhibit selective signatures in the BMX pig genome. In 
addition, future analysis of Asian and European pig genome 
data will help clarify breed development in pigs. 
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ABSTRACT 


Three-finger toxins (TFTs) are well-recognized non- 
enzymatic venom proteins found in snakes. However, 
although TFIs exhibit accelerated evolution, the 
drivers of this evolution remain poorly understood. 
The structural complexes between long-chain 
a-neurotoxins, a subfamily of TFTs, and their nicotinic 
acetylcholine receptor targets have been determined 
in previous research, providing an opportunity to 
address such questions. In the current study, we 
observed several previously identified positively 
selected sites (PSSs) and the highly variable 
C-terminal loop of these toxins at the toxin/receptor 
interface. Of interest, analysis of the molecular 
adaptation of the toxin-recognition regions in the 
corresponding receptors provided no statistical 
evidence for positive selection. However, these 
regions accumulated abundant amino acid variations 
in the receptors from the prey of snakes, suggesting 
that accelerated substitution of TFTs could be a 
consequence of adaptation to these variations. To the 
best of our knowledge, this atypical evolution, initially 
discovered in scorpions, is reported in snake toxins for 
the first time and may be applicable for the evolution 
of toxins from other venomous animals. 


Keywords: Three-finger toxins; Nicotinic acetylcholine 
receptor; Driver 


INTRODUCTION 


Three-finger toxins (TFTs) are among the most abundantly 
secreted and effective components of snake venom. They are 
encoded by a large multigene family and show a diversity of 
functional activities (Fry et al., 2003; Kini, 2011). Members in 
this family contain 57—82 residues and four conserved disulfide 
bridges (Endo & Tamiya, 1987), and are identified by three 
loops extending from a central core resembling their namesake 
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three fingers (Fry et al., 2003). 

The long-term evolutionary process that has resulted in the 
high diversity of TFTs conforms to the birth-and-death model 
of multigene family evolution (Nei et al., 1997). This model of 
evolution produces a series of toxins, allowing snakes to adapt 
to a variety of prey and predators (Nei et al., 1997). Snakes 
employ their venom as a weapon to disable prey (primary 
function) and as a defensive tool against predators (secondary 
function) (Kang et al., 2011). Snake toxins and prey are 
likely involved in a co-evolutionary arms race, whereby evolving 
toxin resistance in prey species and novel toxin evolution in 
snake species exert mutual selective effects (Arbuckle et al., 
2017; Calvete, 2017; Jansa & Voss, 2011; Voss & Jansa, 
2012). However, despite evidence of accelerated evolution 
in the TFT family (Sunagar et al., 2013), the cause of this 
evolution remains vague. 

The TFT family members can target various receptors and 
ion channels with high affinity and specificity (Doley et al., 
2009; Kini, 2011). The a-neurotoxins (a-NTXs) of TFTs 
can interact with nicotinic acetylcholine receptors (nAChRs), 
inhibiting postsynaptic membrane ion flow and leading to 
flaccid paralysis (Barber et al., 2013; Chiappinelli et al., 1996). 
Based on the chain length and number of disulphide bridges, 
a-NTXs are usually divided into short or long a-NTXs. Short 
a-NTXs contain 60-62 amino acid residues with four disulfide 
bonds, whereas long a-NTXs contain 66-75 residues and five 
disulfide bonds (Dajas-Bailador et al., 1998). 

nAChRs are pentameric transmembrane proteins of 
ligand-gated ion channels and are formed by different 
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combinations of five subunits, that is, a, B, y, 6 and e (Corringer 
et al., 2000; Karlin, 2002). Each subunit is composed of a large 
N-terminal extracellular domain, which serves as the major 
binding site for the toxins, followed by four transmembrane 
helices and a small C-terminal extracellular domain (Dellisanti 
et al., 2007; Karlin, 2002). nAChRs can be further divided into 
muscular or neuronal types (Corringer et al., 2000; Wang et 
al., 2002). a-NTXs can potently antagonize the a1 subunit of 
heteropentameric muscle nAChRs ((a@1)28176 in fetal muscle 
and (01),B1e0 in adult muscle) and the a7, a8, a9 or ato 
subunits of homopentameric neuronal nAChRs (Karlin, 2002; 
Sine et al., 2013). Crystal structures in the complexes between 
long a-NTXs and related receptors have been identified and 
offer opportunities to explore the molecular mechanism driving 
the accelerated evolution of these toxins (Bourne et al., 2005; 
Dellisanti et al., 2007; Huang et al., 2013). 

In this work, we found several sites previously identified as 
positively selected sites (PSSs) as well as the highly variable 
C-terminal loop of the toxins to be located at the toxin/receptor 


binding interface (Bourne et al., 2005; Sunagar et al., 2013). 


Structural and evolutionary analyses of the toxin-recognition 
region from the receptors (a1, a7, a9 and «10 subunits from 
nAChRs) uncovered an atypical evolution between snakes and 
their prey, in which the amino acid diversity of the nAChR 
toxin-binding regions appeared to drive the adaptive evolution 
of the TFT family. This study showed good agreement with 
our previous research on scorpion toxins and sodium channels 
from their competitors (Zhang et al., 2015). Furthermore, this 
paper provides a broader vision into the evolution between 
venomous animals and their prey/predators. 


MATERIALS AND METHODS 


Sequence analysis 

ClustalX (http://www.clustal.org) was used to align all 
nucleotide and amino acid sequences. A total of 130 long 
a-NTX sequences were aligned and used for sequence logo 
analysis with the WebLogo program (Supplementary Figure 
S1). For the receptors, 76 sequences of muscle-type a1 
nAChRs (three from Reptilia, three from Amphibian, four from 
Aves, 31 from small mammals, and 35 from fishes), 59 
neuronal-type a7 nAChRs (one from Gastropoda, two from 
Arthropoda, three from Reptilia, three from Amphibian, three 
from Aves, 18 from small mammals, and 29 from fishes), 
68 neuronal-type a9 nAChRs (three from Reptilia, three from 
Amphibian, four from Aves, 30 from small mammals, and 28 
from fishes), and 52 neuronal-type a10 nAChRs (two from 
Reptilia, two from Amphibian, four from Aves, 24 from small 
mammals, and 20 from fishes) (Supplementary Figures S2-S5) 
were also aligned. For positive selection analyses, aligned 
nucleotide sequences of the receptors were used to construct 
neighbor-joining trees. 


WebLogo and ConSurf analyses 


WebLogo can be used to generate sequence logos, which 
are graphical representations of the patterns within multiple 
sequence alignments. The alignments of the aforementioned 
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toxins were used to create sequence logos to identify the 
conservation of each position (Supplementary Figure $1). 
Each logo comprises stacks of letters, one stack for each 
position in the sequence (Crooks et al., 2004). The overall 
height of each stack shows the sequence conservation of that 
position, whereas the height of the symbols within the stack 
indicates the comparative frequency of the amino acid at the 
position (Crooks et al., 2004). Generally, sequence logos 
provide a richer and more precise description of conserved 
and variable regions within sequences. We used the ConSurf 
program (http://consurf.tau.ac.il/) under default parameters to 
calculate conservation scores of the amino acid sequences of 
the related receptors. ConSurf not only depends on sequence 
alignments but also on phylogenetic trees to identify conserved 
and variable regions (Landau et al., 2005). A Bayesian tree 
was constructed using the corresponding alignments with the 
JTT evolutionary substitution model. One of the advantages of 
ConSurf compared with other methods is that the computation 
of the evolutionary rate is more precise when employing the 
empirical Bayesian or maximum-likelinood methods. 


Positive selection analysis 

Excess nonsynonymous substitutions compared with 
synonymous substitutions (w=d\/ds>1) is an important 
sign of positive selection at the molecular level (Yang, 1998; 
Zhu & Gao, 2016). To perform such analysis, we compared 
two pairs of site models (Mia (neutral)/M2a (selection) and 
M7 (beta)/M8 (beta & w)) to measure the selective pressure 
of the receptors to which the long a-NTXs bind (a1, a7, a9 
and «10 subunits from nAChRs). Model M2a and M8 add a 
site class to Mia and M7, respectively, with the free œ ratio 
calculated from the data and used to determine the probability 
of positive selection (Anisimova et al., 2001; Anisimova et al., 
2003; Wong et al., 2004; Yang € Swanson, 2002). As Mia 
and M7 are nested within their respective alternative models 
(M2a and M8) and have two more parameters, x? distribution 
can be used for the likelinood ratio test to compare the fit 
of the two competing models. We used the Bayes Empirical 
Bayes method to calculate the posterior possibility that each 
codon is from the site class of positive selection. The Bayes 
Empirical Bayes method is an improvement of the previous 
Naive Empirical Bayes method and accounts for sampling 
errors in the maximum-likelihood estimates of parameters in 
the model (Nielsen & Yang, 1998). Sites with a high possibility 
(>95%) of coming from the class with œ>1 are likely under 
positive selection and can be analyzed further (Yang, 1998). 


RESULTS 


PSSs of toxins locating on the toxin-receptor complex 
interface 

The N-terminal extracellular domains of nAChRs, which consist 
of a 10 stranded B-sandwich and an N-terminal a-helix, act 
as the major binding sites for long a-NTXs (Changeux et al., 
1970). The three regions of long a-NTXs comprising fingers | 
and Il and the C-terminus are involved in interactions with the 
receptors, with finger Il being the main stabilizing interaction 


center (Bourne et al., 2005). The a211 structure (mouse 
nAChR a1 subunit (PDB entry 2qc1)) can be seen as a 
representative of the N-terminal extracellular domain of the 
nAChR subunit (Dellisanti et al., 2007), in which loops B4-B5 
(loop A), B7-B8 (loop B) and B9-B10 (loop C) serve as the 
principal ligand-binding interfaces (Brejc et al., 2001; Unwin, 
2005) and loop C is the most important region for high affinity 
with long a-NTXs, as revealed by site-directed mutations 
(Fruchart-Gaillard et al., 2002; Levandoski et al., 1999). Loop C 
of «211 is enveloped by fingers | and Il and the C-terminal loop 
of the toxin, with finger ll inserted into the ligand-binding site 
wrapped by loops A, B and C of a211 (Figure 1A) (Dellisanti 
et al., 2007). In addition, finger | is sandwiched by loop C, 
whereas finger Ill weakly contributes to «211 binding (Dellisanti 
et al., 2007). 
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Figure 1 a-bungarotoxin interacts with «211 and the sequence 
logos of long a-NTXs 

A: Interaction between «a-bungarotoxin (long a-NTXs) and a211. 
The PSSs 


of the long a-NTXs are in pink. The three loops of #211, major toxin-binding 


a-bungarotoxin is shown in wheat and «211 in light blue. 
regions, are highlighted in blue. B: Sequence logos of long a-NTXs. The 
C-terminal loop of the long a-NTXs involved in receptor binding is indicated 


in a red box. 


We mapped the PSSs of long a-NTXs identified by Sunagar 


et al. (2013) on the complex structure of a-bungarotoxin and 
its receptor (Figure 1A). The toxin-receptor complex model 
showed that some PSSs of the long a-NTXs were located 
on the toxin binding surface with the receptors (Bourne et al., 
2005; Dellisanti et al., 2007; Huang et al., 2013). Almost 
all PSSs in finger Il (Ala31, Phe32, Ser34, Ser35 and 
Val39) of a-bungarotoxin are involved in the interactions with 
the receptors (Dellisanti et al., 2007; Dimitropoulos et al., 
2011; Huang et al., 2013). The same situation appears in 
a-cobratoxin (long a-NTX) since their fast evolved sites (Ala28, 
Phe29, Ser31, lle32 and Arg36) in finger Il also overlap with 
the toxin binding sites of a-cobratoxins (Bourne et al., 2005; 
Dimitropoulos et al., 2011). 

The C-terminal loop of the long a-NTXs also contributes 
to the interaction with related receptors. Multiple sequence 
alignments of long a-NTXs, generated by ClustalX, were used 
to create sequence logos (Figure 1B, Supplementary Figure 
S1) (Crooks et al., 2004). The C-terminal loop of the long 
a-NTXs, indicated in the red box in Figure 1B, was highly 
variable based on the WebLogo analyses, and the whole C-tail 
involves extensive insertions and deletions. Thus, this loop 
might undergo positive selection despite the technical difficulty 
in detecting PSSs from indel-containing sequences. 


No evidence for positive selection in the toxin-recognition 
regions of receptors 


Snakes employ their venom to immobilize various prey, 
including snails, insects, fishes, toads, lizards, chickens, small 
mammals and even other snakes (Kang et al., 2011). The 
maximum-likelihnood models of codon substitutions were used 
to identify selective pressure in the toxin-recognition regions 
of the receptors from the prey of snakes. However, no 
positive selection signals were detected in the receptors of long 
a-NTXs (a1, «7, 019 and 0,10 subunits) (Table 1, Supplementary 
Tables S1—S4). The maximum-likelinood estimates under MO 
showed that the average œ ratios for all receptor sequence 
pairs ranged from 0.02 to 0.07. M2a and M8 detected no 
evolution-accelerating sites. The ws under M8 of the a1 and 
&œ10 subunits of nAChRs from snake prey equaled 1 (Table 
1). Under the M8 model, the «7 and 0.9 subunits of nAChRs 
showed (ws>1 (Table 1), but their proportions (p1) equaled O, 
proving that no PSSs existed (Supplementary Tables 81-84). 


Table 1 Parameter estimates and likelihood ratio statistics (2A1) for different subunit types of nAChRs 


Type l (M1a) l (M2a) wa (M2a) 
a1 subunit  —12874.0 -—12874.0 1.00 
a7 subunit  —13163.9 —13163.9 1.00 
a9 subunit  —11856.3 -11856.3 1.00 
0110 subunit —9954.2 —9954.2 1.00 


2Al 1(M7) 1(M8) ws (M8) 2A 


O -12617.4 -—12607.4 1.00 20 
O -12785.7 -—12785.7 34.66 0 
O  —11626.9 -—11626.9 1.42 0 
0 -9746.8  -—9746.8 1.00 0 


l is the log likelihood; 2A/ is between null models and their alternative models: M1a/M2a and M7/M8. Two o values in M2a and M8 are 


boldfaced. 
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Snake toxins bind to variable regions of their receptors 
Although PSSs were detected in the toxins, none were found in 
the toxin-recognition regions of the involved receptors. Based 
on the complex models between the long a-NTXs and related 
receptors, we further analyzed the evolutionary conservation 
of the N-terminal extracellular domain regions of the nAChR «a1 
and a7 subunits from snake prey using ConSurf (Dellisanti et 
al., 2007; Huang et al., 2013) (Figure 2). Our results indicated 
that loop C demonstrated the greatest variation among the 
three receptor loops (Figure 2). 
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Figure 2 ConSurf results of «-bungarotoxin with the receptors 
ConSurf results of a-bungarotoxin with «1 (PDB entry 2qc1) and a7 (PDB 
entry 4hqp) nAChRs are shown in A and B, respectively. ConSurf provides 
evolutionary conservation scores for nAChRs. PSSs of finger | and finger II 
and the residues of the C-terminal loop in a-bungarotoxin are highlighted in 
red (PSSs of finger II are partly hidden). Toxin-binding regions (loop A, loop 
B and loop C) of the a1 and a7 nAChRs are spherized. Binding sites variable 


in loop C are indicated. 


434 www.zoores.ac.cn 


We analyzed the interaction pairs in the «a-bungarotoxin 
and a1 nAChR subunit complex, which were lle11-Phe189/ 
Pro194, Val31-Tyr93/Asp99/Phe100, Phe32-Tyr93/Phe100/ 
Trp149/Tyr190, Val39-Val188/Tyr190, His68-Pro194, Pro69- 
Ser191, Lys70-Pro194 and Gln71-Pro194, with the toxin sites 
corresponding to the PSSs in fingers | and Il and the C-terminal 
loop (Dellisanti et al., 2007; Dimitropoulos et al., 2011). The 
loop A (Tyr93, Asp99 and Phe100), loop B (Trp149) and loop 
C (Val188, Phe189, Tyr190, Ser191, Pro194 and Pro197) 
were involved in the interaction with the PSSs of the toxin 
binding surface. Although Tyr93 in loop A and Trp149 in loop B 
were conserved, Asp99 and Phe100 in loop A exhibited some 
variability (99: Asp/Asn and 100: Phe/Tyr). In contrast, the 
binding sites in loop C were more variable (188: Val/Arg/Lys, 
189: Phe/Tyr/Thr, 191: Ala/Ser/Thr/Pro and 194: Pro/Gln) 
(Figure 2A). 

For the a-bungarotoxin and a7 nAChR subunit complex, 
the interaction pairs include lle11-Phe183/Lys188, Phe32- 
Tyr91/Trp145/Arg182/Tyr184, Val39-Arg182/Tyr184, His68- 
Lys188, Pro69-Glu185, Lys70-Cys186/Lys188 and Gln71- 
Lys188, with the sites of the toxins also corresponding to 
the PSSs in fingers | and Il and the C-terminal loop (Huang 
et al., 2013). Compared with the conserved residues in loop 
A and loop B (Tyr91 and Trp145), the binding sites in loop 
C were diverse (182: Lys/Arg/Leu/Ser/Asn, 183: Phe/Tyr/lle, 
185: Glu/Asp/Asn/Gly and 188: Lys/Asp/Glu/Pro/Gin) (Figure 
2B). The additional PSSs (Val31, Ser34 and Ser35) in finger 
Il contacted residues of the complementary subunit, and 
included Ser32, Ser34, Leu36, Trp53, Gln55, Gln114, Leu116 
and Asp160 (Huang et al., 2013). These residues in the 
complementary subunit may provide the driving force for the 
additional sites in finger Il. 

Taken together, our results suggest that the variable 
toxin-recognition region in the receptors might drive the 
accelerated evolution of the toxin-binding residues. 


DISCUSSION 


By mapping the PSSs of the long a-NTXs on the toxin-receptor 
complex, we found that several of them are located on the 
toxin-binding surface of the receptor (Figure 1A). In addition 
to these PSSs, high sequence diversity was also observed 
in the C-terminal loop of the long a-NTXs. Thus, given its 
importance in the interaction with receptors, we surmised that 
it could be an accelerated substituted region for adaptation to 
receptor variability (Figure 1B). Several PSSs were detected 
in finger IIl of the long a-NTXs, although finger III contributed 
little to the binding. This may be due to its high flexibility 
in the complex. Other mechanisms may be involved in the 
accelerated substitutions of this finger, which requires further 
investigation. 

We further observed the amino acid variability of the 
principal toxin-recognition regions (mainly in loop C) in related 
nAChR subunits (Figure 2). Compared with loop A and loop B, 
loop C exhibited the greatest variation due to its predominant 
role in the toxin-receptor interactions. Thus, we concluded that 
the evolutionary variability of the toxin-recognition regions of 


the receptors is a possible driver for the accelerated evolution 
of the toxins. 

Short-chain a-NTXs can bind to nAChRs with high affinity 
(Trémeau et al., 1995). Loop II of short-chain a-NTXs pushes 
into the ligand-binding pocket of nAChRs, whereas the tips 
of loops | and Ill contact nAChRs only in a ‘surface-touch’ 
way (Mordvintsev et al., 2005). Previous study on Torpedo 
californica showed that the C-loop is vital for the binding of 
short a-NTXs to nAChRs (Bourne et al., 2005; Mordvintsev et 
al., 2005). However, as there are no specific coordinate files of 
complexes between short a-NTXs and their targets, hindering 
further study. 

Our previous study on the evolution of scorpion toxins and 
voltage-gated sodium (Nay) channels indicated that variability 
of the toxin-recognition regions in the Nay channels from 
scorpion predators and prey is a putative driver of the 
accelerated evolution of the functional regions of scorpion 
toxins following gene duplication (Zhang et al., 2015). Similarly, 
PSSs exist in the binding interface of the long a-NTXs, 
though only amino acid variability was detected in the 
principal toxin-recognition regions of the related nAChR 
subunits, suggesting that amino-acid substitutions on the 
toxin-recognition surface in related nAChR subunits could 


provide a force driving the accelerated evolution of the toxins. 


Thus, the atypical co-evolutionary manner between snake 
toxins and their receptors is similar to our previous research 
on scorpion toxins and their targets (Zhang et al., 2015). We 
supposed that accelerated evolution of the receptor-bound 
regions of the snake toxins is a consequence of adaptation 
to variable receptors of their prey. From the viewpoint of 
a co-evolutionary arms race between predators and prey 
(Arbuckle et al., 2017; Calvete, 2017; Jansa & Voss, 2011; 
Voss & Jansa, 2012), it appears that prey might exert greater 
selective pressure on their predators, as described in our 
current study. As this evolutionary manner has been shown 
to occur in two distant species, we believe it will be revealed in 
more toxins from diverse venomous animals. 
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